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I.  Summary 


The  research  performed  under  the  contract,  during  the  period  14  July  1986  through 
13  January  1987,  can  be  divided  into  two  main  topics;  coupling  line  source  calculations 
in  2-D  lateral  variations  of  near  source  structure  to  analytic  teleseismic  body  wave  codes 
,  and  the  effects  of  ocean  continent  transition  zones  on  Lt  waves. 

In  section  II,  we  present  a  new  method  for  interfacing  numerical  and  integral  tech¬ 
niques  which  allows  greater  flexibility  in  seismic  modeling.  Specifically,  numerical  calcula¬ 
tions  in  laterally  varying  structure  are  interfaced  with  analytic  methods  that  enable  pro¬ 
pagation  to  great  distances.  Such  modeling  is  important  in  studying  situations  contain¬ 
ing  localized  complex  regions  not  easily  handled  by  analytic  means.  The  calculations 
involved  are  entirely  two-dimensional,  but  use  of  an  appropriate  source  and  then  a  filter 
applied  to  the  resulting  seismograms  produce  synthetic  seismograms  which  are  point 
source  responses  in  three  dimensions.  The  integral  technique  is  called  two-dimensional 
Kirchoff.  Data  from  Yucca  Flat,  Nevada  Test  Site  (NTS)  are  modeled  as  a  demonstra¬ 
tion  of  the  usefulness  of  the  new  method.  In  this  application,  both  local  and  teleseismic 
records  records  are  modeled  simultaneously  from  the  same  model  with  the  same  finite- 
difference  run.  The  application  documents  the  importance  of  locally  scattered  Rayleigh 
waves  in  the  production  of  teleseismic  body  wave  complexity  and  coda. 

In  section  III,  the  effects  of  the  length  of  the  intermediate  path  between  the  con¬ 
tinent  ocean  and  ocean  continent  transition  regions  is  investigated.  First  the  results  of 
two  FE  calculations  with  different  intermediate  path  lengths  are  presented  and  com¬ 
pared.  These  examples  are  contrasted  with  the  path  length  used  in  the  previous  report. 
Then  the  RT  integration  method  is  discussed  and  explained.  First  analytic  expressions 
for  the  stress  components  of  double  couple  and  line  sources  are  derived,  then  the  expres¬ 
sions  for  displacement  and  stress  line  source  Green’s  functions  are  determined.  These 
expressions  are  used  to  illustrate  the  validity  and  determine  the  accuracy  of  the  RT  cou¬ 
pling  method.  The  RT  coupling  method  is  then  used  to  continue  the  propagation  of  FE 
results  through  a  layered  structure  using  the  displacement  and  stress  Green’s  functions 
for  the  remaining  path  length  and  the  displacement  and  stress  time  histories  recorded  at 
a  column  of  element  centers  within  the  FE  grid.  The  total  attenuation  of  Lt  due  to  pro¬ 
pagation  through  a  2-D  transition  of  continent  to  ocean  to  continent  is  found  to  be  at 
most  a  factor  of  four  to  six.  This  is  inadequate  to  explained  the  observed  attenuation  of 
Lt .  Thus,  additional  effects,  other  than  geometry  must  be  considered  to  provide  a  com¬ 
plete  explanation  of  the  attenuation  of  Lt 


Numerical- analytical  interfacing  in  two  dimensions 
with  applications  to  modeling  NTS  seismograms 


Richard  J.  Stead  and  Donald  V.  Htlmhtrger 
Seismological  Laboratory 
California  Institute  of  Technology 
Pasadena,  CA  91125 


March  27,  1986 


submitted  to  ” Scattering  and  Attenuation  of  Seismic  Waves ”, 
special  issue  of  Pure  and  Applied  Geophysics 


I 


suggested  short  title: 

2d  numerical-analytical  interfacing 


■  -  r  t  *.v 


v 


\t ->  v  w\, >ru rv 


-2- 


ABSTRACT 

A  new  method  for  interfacing  numerical  and  integral  techniques  allows 
greater  flexibility  in  seismic  modeling.  Specifically,  numerical  calculations  in 
laterally  varying  structure  are  interfaced  with  analytic  methods  that  enable  pro¬ 
pagation  to  great  distances.  Such  modeling  is  important  in  studying  situations 
containing  localized  complex  regions  not  easily  handled  by  analytic  means.  The 
calculations  involved  are  entirely  two-dimensional,  but  use  of  a  appropriate 
source  and  then  a  filter  applied  to  the  resulting  seismograms  produce  synthetic 
seismograms  which  are  point  source  responses  in  three  dimensions.  The  integral 
technique  is  called  two-dimensional  Kirchhoff  because  its  form  is  similar  to  classic 
three-dimensional  Kirchhoff.  Data  from  Yucca  .Flat,  Nevada  Test  Site  are 
modeled  as  a  demonstration  of  the  usefulness  of  the  new  method.  In  this  applica¬ 
tion,  both  local  and  teleseismic  records  are  modeled  simultaneously  from  the 
same  model  with  the  same  finite-difference  run.  This  application  documents  the 
importance  of  locally  scattered  Rayleigh  waves  in  the  production  of  teleseismic 
bodywave  complexity  and  coda. 


Introduction 

Seismologists  have  long  recognized  that  structural  complexities  near  seismic 
sources  may  affect  teleseismic  waveforms.  For  instance,  events  occurring  in  the 
Imperial  Valley,  California  produce  extended  teleseismic  short  period  signals  last¬ 
ing  much  longer  than  expected  from  near-in  strong  motion  studies,  see  Hartzell 
and  Heaton  [1982].  Presumably,  the  energy  trapped  by  the  low-velocity  layering 
scatters  out  the  bottom  of  the  basin  when  it  encounters  the  basin  edge. 
Recently,  Vidale  and  Helmberger  [1987a],  using  a  finite-difference  scheme,  had 
considerable  success  at  modeling  a  profile  of  the  San  Fernando  earthquake  strong 
motion  records  that  cross  the  Los  Angeles  basin.  Their  approach  assumes  two- 
dimensional  symmetry,  as  displayed  in  Figure,  1,  but  corrects  for  three- 
dimensional  spreading  and  mimics  the  well  known  double  couple  radiation  field. 
Essentially,  the  numerical  excitation  is  matched  to  an  asymptotic  analytical 
source  representation.  Synthetics  generated  by  this  procedure  match  closely 
those  generated  by  analytical  point  source  codes  for  the  same  flat  layered  case 
[Vidale,  et  al.,  1985].  However,  these  solutions  cannot  be  propagated  to  great  dis¬ 
tance  because  of  computational  efforts.  Thus,  a  technique  is  needed  to  interface 
the  numerical  field  back  into  an  analytical  scheme  such  that  the  signals  can  be 
sent  to  large  distances.  This  is  the  basic  objective  of  this  paper,  and  we  will  also 
discuss,  as  a  demonstration,  a  well  controlled  experiment  where  scattered  surface 
waves  can  be  seen  leaving  the  local  field  to  reappear  as  teleseismic  body  waves. 

The  most  controlled  experiment  consistent  with  the  above  motivation  above 
is  that  of  explosions  fired  at  Yucca  Flat  in  the  Nevada  Test  Site.  The  local  struc¬ 
ture  of  Yucca  Flat  is  a  basin  containing  volcanic  tuffs  and  alluvium  [Eckren,  1968 
and  Houser,  1988].  The  events  in  this  area  have  a  striking  complexity  on  telese¬ 
ismic  records  for  an  explosion  source  (see  Figure  2).  Various  studies  show  the 
structure  of  this  region  to  be  seismically  complex  [Taylor,  1983,  among  many  oth¬ 
ers].  Hart,  et  al.  [1979]  examine  the  variation  of  amplitudes  and  magnitudes 
within  Yucca  Flat  dependent  on  source  position  in  the  valley.  Recent  studies 
have  concentrated  on  the  azimuthal  variations  observed  in  the  data  in  both  the 


time  [Lay,  Wallace  and  Helmberger,  1984]  and  frequency  [Lay,  1986a]  domains. 
Studies  of  these  records  have  indicated  the  presence  of  local  scattering  structure 
[Lay,  1986b].  In  addition,  existing  strong  motion  records  demonstrate  lateral 
anomalies  in  the  propagation  of  seismic  energy  at  Yucca  Flat.  Figure  3  shows 
the  vertical  velocity  records  from  the  event  FLASK.  The  large  difference  in 
amplitude  and  duration  of  both  the  first  arrivals  and  the  Rayleigh  waves  from 
east  to  west  strongly  suggests  that  scattering  plays  an  important  role.  It  is  found 
herein  that  these  scattered  Rayleigh  waves  are  the  likely  progenitors  of  telese- 
ismic  complexity  apparent  in  Figure  2.  This  is  consistent  with  the  results  of 
Lynnes  and  Lay  [1987],  published  in  this  issue. 

A  Two-Dimensional  Representation  Theorem 

In  accordance  with  the  motivation  discussed  above,  we  derive  a  two- 
dimensional  representation  theorem  method  similar  in  its  application  to  three- 
dimensional  Kirchhoff  integration.  Such  an  integral  method  allows  computation¬ 
ally  inexpensive  generation  of  teleseismic  records  from  complex  source  region  cal¬ 
culations,  particularly  those  of  finite  difference.  In  deriving  such  a  method,  the 
resulting  expression  must  be  fully  elastic,  require  a  minimum  of  computation 
time,  and  be  readily  adaptable  to  finite-difference  methods.  A  straightforward 
approach  is  to  parallel  the  derivation  of  three-dimensional  integration  methods. 
For  more  information  on  such  integral  methods  in  both  two  and  three  dimensions 
see  Baker  and  Copson  [1950]  and  Mow  and  Pao  [1971]. 

Starting  with  the  two-dimensional  elastic  wave  equation  in  polar  coordinates, 

d3u  _I_  _l_  1  /.  % 

dr 3  r  dr  +  r*  df  ”  aa  dt 3 

take  the  Laplace  transform  over  time,  ignoring  the  6  -  dependence,  and  setting 


1  " -F  -"Ul  'VW  -  tf  *  V  *  V*  m  IT*  IfU  L~V  VVtTVl  ^W*’\**VT/"mV  7/  rtw .V^  kTW  VHk 


which  has  solutions  like  /f0|  — j,  see  Hudson  [1963]  for  example. 


Thus,  u  depends  only  on  r  and  solves 
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Let  w  be  another  solution  of  the  Equation  2,  then  Green’s  transformation 


becomes 
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The  geometry  for  this  transformation  is  shown  in  Figure  4.  Now  let 
w  Then  for  P  £  D, 

For  P  6  D,  as  in  Figure  5,  the  log  singularity  at  P  for  gives  the  result 
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Let  n  =  -~/r  in  the  integral  over  <r,  then 
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where  dc  is  a  length  element  on  the  circumference  of  the  inner  circle.  Again,  let 
w  =  K o,  and  also  let  the  radius,  t ,  of  a  go  to  zero. 
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Near  P,  K0  can  be  approximated  as  follows. 


K0  — I  ~  -log 

O  )  2  (  Or 


(8) 

(9) 


If  these  approximations  are  substituted  into  the  Equation  7,  the  integral  over  a  is 
evaluated  in  the  limit. 
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Therefore, 
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This  integral  can  be  exactly  inverted  to  the  time  domain.  First,  take  the 


derivative  of  /Co- 


Substitute  this  into  the  integral. 


Take  the  inverse  transform. 
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Inverting  the  above  integral  requires  the  use  of  the  convolution  rule,  the 


derivative  rule  and  the  inverse  transforms 


L-\K<lbp ))  [6  >0]  -  {tUy*  H{t-b) 
L-\KAbp))  [»  >0]  -  |  (tUy*H(t-b ). 
Now,  taking  the  inverse, 

For  simplicity,  let  /(i )  —  (<Mr  /£»)*)‘1/,^(Mr  /<*)• 


dn  [  Vt'-ir/af 


•|p,~^r/r[7l;‘^‘(,'/(‘»+lr'/(')] '  <17> 

Equation  17  may  be  cast  in  a  form  similar  to  that  of  conventional  Kirchhoff 
methods.  First  apply  the  convolution  identity  /(f)*? (t  +e )*“/(t+c )*g (t)  with 


c  —  -r  /a. 


«(P)-  JLf 

1  '  2ir  Jr  [r  an  it  U  a 


)/(<))  +  fr/(<)  41 


where  «  is  now  «(r),  r—  f-r/a,  the  retarded  time,  and  the  operator 

/(l>"  77  "+\rla-  R'WrUm*  th“’ 

Bring  the  derivative  across  the  convolution  in  the  last  term. 

-  H  ['“>•(£■ ♦  ifcSHwrvfSSk +  *ffc]  - 

Rearranging  terms  and  taking  the  last  time  derivative, 

*ss  * 


I 

< 


This  equation  is  an  exact  form  of  the  representation  theorem. 


It  is  useful  to  consider  an  approximation  of  Equation  18  for  large  r.  The 
main  reason  for  considering  this  is  that  /(<)  depends  on  r  and  therefore  the  con¬ 
volution,  which  is  computationally  expensive,  cannot  be  moved  outside  the 
integral.  For  large  r,  the  operator  f(t)  is  dominated  by  l/\/T ,  that  is 

/(()  =  - - -  r ,and 

'  '  yjt  y/t+lr/a  vtv2rja 

T-±-7-f(t)~H(t)VT(2r/a)>* 

Equation  18  then  becomes 


«(P) 


1  sfot  1  f  1  f  «  dr  )  JI 

2 it  y/i  vT  Jr  v/r"  l  dn  2 r  3n  a  dn  dt  J 


2jt  \/2  Jr  2r  Vr  2r  3n 

The  second  integral  is  similar  to  the  second  term  in  the  first  integral,  except  that 
it  falls  off  as  l/r  with  respect  to  it.  This  means  that  the  second  integral  may  be 
ignored  at  large  r . 


_1_  y/o  1  t  j"  If  du_  «  dr  1  dr  dn  | 

2ir  v/2  v/T  <J r  l  9n  2r  dn  a  dn  dt  J 


<0 


(19) 


This  approximation  also  arises  from  using  the  asymptotic  form  of  K0  before 
inverting  the  transform.  That  is, 


K0 


(2P) 


Substitute  Equation  20  into  Equation  11,  taking  the  appropriate  derivatives,  to 


get 


(21) 


This  may  be  inverted  using  the  shift  rule  to  account  for  the  <  * ,  and 
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L~l(\/irJT )  —  1/vT,  to  arrive  at 


«(P)  — 


1  1  1  t  I"  £  j  dv_  «  3r  1  dr  dv  \ 

2jt  \/2  n/T  Jr  v/7  (  3n  2r  3n  cr  3n  3(  J 


dl. 


Once  again,  «  is  a  function  of  retarded  time  r  =  t-r/a,  this  time  due  to  the 
application  of  the  shift  rule. 

It  is  interesting  to  compare  Equation  19  to  the  more  familiar  point  source 
Kirchhoff  formula  from  optics  [see,  for  example,  Baker  and  Copson,  1950). 


«(P)~ 


1  dr  du  1  3a  | 
ar  dn  dt  r  dn  ) 


which  can  be  written 


«{P) 


1  I"  r  -1_  3a  ^  dr  1  dr  3a 
4 n  J  Js  r  dn  r  dn  a  dn  dt 


dS 


(22) 


(23) 


Comparing  Equations  19  and  23,  the  form  is  identical  between  the  integrals, 
except  for  a  factor  of  two  on  the  second  term.  The  convolution  with  1/vT  and 
the  difference  in  scaling  with  distance  are  expected  for  a  line  source  as  opposed  to 
a  point  source.  Both  integrals  can  be  applied  in  similar  ways.  Given  a  two- 
dimensional  problem,  the  method  of  choice  should  be  the  line  source  integral. 


An  analytic  verification  of  these  formulae  in  a  whole  space  demonstrates  the 
operation  of  both  the  three-dimensional  Kirchhoff  formula  and  the  two- 
dimensional  integral  just  derived.  The  verification  of  the  three-dimensional  Kir¬ 
chhoff  formula  is  similar  to  that  of  Hilterman  [1975]  who  set  up  a  reflected  case. 
The  setup  of  the  problem  is  given  in  Figure  6,  where  a  plane  (S)  is  assumed 
located  equidistant  between  the  source  and  receiver  (that  is,  r0=r).  Conceptu¬ 
ally,  the  source  lights  up  the  surface,  and  the  surface  reradiates  the  energy  to  the 
receiver.  In  this  case,  the  transmission  coefficient  is  unity  so  that  only  geometri¬ 
cal  effects  are  tested.  Starting  with  the  integral  in  Laplace  space, 


ftp)  = 


_L  ff  Jr 

1  3<£(S)  <t>{ S)  dr  t  d>(S)  dr 

4»r  J  J s 

r  dn  r2  dn  a  r  dn 

(24) 


In  this  example,  #p)  is  the  potential  at  the  observation  point,  and  #s)  is  potential 
on  the  interface  S .  For  a  symmetric  point  source, 


— «  *  /(») 
ro 


Applying  the  chain  rule  to  ~s~, 

on 


asa  .  «sa  *•_.  i  i+j.  »*•*,,  (28) 

dn  3r0  dn  {  a  r0  )  dn  v  ' 

Substituting  Equations  25  and  26  into  Equation  24  yields 

*» -  *•  *JJ. (*♦*) w 

Let  r0— r,  4^.— >-4^-  and  ■  where  0  is  solid  angle.  Then 

dn  on  r*  on 


Let  t  =2r /a,  to*,2ro/°  and  d(i^--^-dt . 


Invert  to  the  time  domain. 


1 

i  d3n 

_2_ 

2jt 

a  dt3 

o( 

The  solid  angle  swept  out  over  the  surface  as  a  function  of  time  is 


*£  =  2n±H(t-t0) 


iSL  =  2>r-|y<5{*  -*o)  -  (*  ~*o) 

Ol  (  I 
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Using  these  expressions,  Equation  28  becomes 


^(P)“  2^[2V^"‘o)_4V^(<"to)+  T2,r^(t~,o))*/(0 


Or  t 

1  /(t-to) 

a  t  o 


=  l/(t-<0)  (29) 

This  is  the  familiar  geometric  spreading  law  for  a  point  source. 

The  same  analysis  can  be  performed  in  two  dimensions  (see  Figure  7).  The 
approximate  form  of  the  integral,  as  given  in  Equation  19,  is  used  because  the 
form  is  more  like  that  of  the  three-dimensional  case,  and  the  analysis  is  more 
straightforward.  For  an  appropriate  correspondence,  replace  «  inside  the  integral 
of  Equation  19  with  *(rj,  the  potential  on  the  contour  r.  As  above,  starting  with 
the  Laplace  transform  of  the  integral  (Equation  21),  use  the  approximate  form  of 
the  source  (see  Equation  20). 

4>(n  =  ■  t  °  /(«)  (30) 

V“,r0 

ddXH  Vna  (  1  .  *1  dt\ 


Substitute  these  into  Equation  21  to  obtain 


JZ  f  .-r 


Let  rD-r,  and  d  e=——dl ,  where  $  measures  the  angle  from  vertical, 

dn  dn  r  dn 

as  indicated  in  Figure  7.  Then 


'>  >7 


a  2 r  I  dn 
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Let  f=2r/ar,  tQ=2rQ/a  and  d  $=-—dt . 

■» 


m- /<•)}/•-  («+£)•£" 


Invert  to  the  time  domain. 


*(P,~?  4f  +  / ‘777rfr  *(‘-‘o)'/(0 

o 

The  total  angle  swept  out  over  the  contour  as  a  function  of  time  is 

tf(0=2co.-1[il]/f(r-t0) 
d8  _  „‘o  H(t-t o) 

*  "  «  V&TT 

Using  these  expressions,  Equation  33  becomes 


<o  H(t-to)  f  l  *o  #(!-«>)  , 

<  ^ij  J  »  '  7  "(  *)/(  0 


t0  ,  ,  f  i  r 

'  v^7  °{ 


+  hj "3'/T_  Tdr  ff(l~to)*f(t) 


t0  1  ,  \f?-tj 

'  v/^  +  ‘o*r 


7T77 

'o'  v 


H{t-t0)*f(t) 


to'+t'-tj 
t0t  y/t3-to2 


H(t-t  oW(0 


*o 


is 


(34) 


m 


H(t-t  o) 


'/(*) 


This  is  the  familiar  equation  for  line  source  response.  Note  that  in  the  three- 
dimensional  point  source  case,  dd/dt—O  until  <*f0,  and  then  jumps  to  2n/t0  at 
gradually  decreasing  with  time  thereafter.  An  interface  with  structure 
behaves  in  a  more  interesting  manner,  see  Scott  and  Helmberger  [1984].  In  two 
dimensions,  dO/dt  has  a  square  root  singularity  at  t=»t0-  In  the  three-dimensional 
extension  of  this  case,  the  response  at  f— 10  represents  the  integrated  energy 
arriving  at  the  receiver  from  an  infinite  strip  represented  by  the  first  line  element. 
It  is  the  integration  over  this  infinite  strip  in  and  out  of  the  plane  of  Figure  7 
that  creates  the  singularity. 

The  simplest  application  of  two-dimensional  ICirchhoff  is  to  a  medium  of 
constant  velocity  set  in  Cartesian  coordinates,  with  part  of  the  contour,  r,  paral¬ 
lel  to  one  of  the  axes  (Figure  8).  Choosing  a  boundary  paralleling  the  x-axis  in 
an  (x  ,z )  system  and  extending  to  infinity  in  both  directions,  the  contour  r  may  be 
closed  at  r  —  oo  in  the  direction  +oo.  As  r-«,  the  integral  over  this  closure 
of  r  vanishes.  The  point  P  becomes  the  receiver  location,  and  all  sources  are  out¬ 
side  the  contour.  The  vector  7  has  its  origin  at  P  and  is  the  distance  to  r.  From 

9  9 

this  information,  n  =  -i  on  the  straight  boundary,  — -  -•  —  and  -r—  —  .  For 

dn  r  on  dz 


this  case,  the  integral  from  Equation  19  becomes 


«(P) 


1  y/a  1  1  (  x  »  dn  | 

2k  v/2  Vt  J  sTr  1  dz  *2 r3  ra  dt  j 

-OO 


dl. 


(35) 


For  computation,  the  integral  must  be  discretized  to  a  finite  sum.  Assuming  that 
the  contributions  from  the  ends  of  the  boundary  are  small, 


h  \/a  1  \  ^  1  (  du  7  7  d*  \  1 

” (P|  ”  ‘2j  [7? I'aT  +  “ +  “  I ) ' 


(36) 


where  h  is  the  spacing  between  discrete  line  elements.  If  the  coordinates  of  P  are 
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tr,  ip,  and  the  coordinates  of  the  first  line  element  are  *„  then 
*,•  »*p-*i-(i-1)A  and  r,  —  \/z?+x*. 

The  above  geometry  is  developed  specifically  for  application  with  finite- 
difference  techniques.  For  the  acoustic  case,  the  wave  field  «  becomes  pressure 
(or  dilatation,  6),  and  Equation  36  is  directly  applicable.  For  the  elastic  case, 
some  modifications  must  be  made.  Assuming  continuous  material  properties 
across  the  boundary,  the  compressional  and  shear  wave  fields  may  be  separated. 
Separating  these  fields  following  a  fully  elastic  finite-difference  calculation 
requires  the  taking  of  the  divergence  and  curl  of  the  full  wave  displacement  field 
for  the  compressional  and  shear  components,  respectively.  This  can  be  seen  from 
the  potential  form  of  the  vector  displacements,  *  —  V^+Vx&  First,  take  the 
divergence,  V  j,  -  and  take  the  curl,  VXju  =  VxVXjfc.  Take  divergence 

and  curl,  respectively,  a  second  time.  Two  equations  now  exist,  one  in  which  *  is 
identified  as  the  second  spatial  derivative  of  the  compressional  wave  field  with 
velocity  a  as  the  p-wave  velocity,  and  a  second  in  which  «  is  identified  as  the 
second  spatial  derivative  of  the  shear  wave  field  with  velocity  a  as  the  s-wave 
velocity,  0.  Returning  to  the  wave  equation,  the  reason  for  the  second  derivatives 
becomes  evident.  The  wave  equations  for  the  second  derivatives  of  displacement 
are 


Vs, 


1 

o*  dt2 


va* 


I  d2!L 
(?  dt 3 


Identifying  the  wave  field  «  in  Equation  10  as  the  second  derivatives  of  the  vec¬ 
tor  displacements  and  j, ,  the  wave  equations  allow  the  second  spatial  deriva¬ 
tives  of  *,(P)  and  ju(P)  to  be  equated  to  second  time  derivatives,  that  is  accelera¬ 
tions.  Thus,  Equation  19  becomes  two  integrals,  one  for  the  p-wave  field  and 
another  for  the  s-wave  field,  in  each  of  which  the  required  inputs  are  spatial 
derivatives  of  displacement  (calculated  in  the  finite-difference  code),  and  the  out¬ 
puts  are  accelerations. 


1 
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Normally,  when  calculating  a  teleseismic  record  in  this  manner,  elastic  or 
acoustic,  a  point  source  in  three-dimensional  geometry  is  the  response  desired. 
An  appropriate  approximate  conversion,  derived  in  Appendix  A,  is 

r"hu^  (A-28> 

Here,  *  is  the  horizontal  distance  (as  in  the  above  derivation),  and  R  is  the  total 
distance  (v'*2+z2)  from  the  source  to  the  receiver.  This  conversion  assumes 
cylindrical  symmetry  about  a  vertical  axis  through  the  source.  If  it  is  applied 
directly  to  Equation  36,  the  result  is  as  follows. 


v/o  2  d  l  „  i  l  (  d%  t 

,(p)  2jt  TiTtr+JTliTT  7T  Ij  [7^|-a7  +  '^T 


£ 

r*a 


This  equation  can  be  simplified  by  recognizing  that 


-  *6(1  WO 


-WO 


Then 


_  _  *  'fa  \ "'  f  1  f  ,  /  /  d*  l  1 

“(P|  _  UTTTTli  2-i  [TTraT  +  *  V  +  7r="a7)  J 


This  conversion  assumes  the  original  line  source  radiation  pattern  was  appropri¬ 
ately  adjusted  for  conversion  to  point  source,  as  discussed  in  Appendix  A. 

Many  assumptions  are  made  in  the  above  derivation.  While  the 
justifications  given  are  reasonable,  some  simple  tests  of  the  method  against  other 
analytical  methods  is  useful.  Suppose  we  examine  a  simple  half-space  with  a 
point  source  at  a  depth  of  1  km  in  an  elastic  half  space  with  p- velocity  4  km/s, 
s-velocity  2.3  km/s  and  density  2.7  g/cc,  and  receivers  at  a  depth  of  6  km  and  a 
variety  of  ranges.  A  comparison  of  the  new  method,  ^finite-difference  and 


Cagniard  seismograms  is  shown  in  Figure  9.  The  seismograms  shown  include  an 
RDP  source  convolved  with  a  Gaussian.  The  Gaussian  is  necessary  in  the  finite- 
difference  scheme  to  limit  the  bandwidth  of  the  source  to  a  range  in  which  propa¬ 
gation  is  stable.  The  Gaussian  is  normalized  to  have  an  integral  of  unity.  The 
records  shown  are  all  velocity  records,  representative  of  the  response  of  a  broad 
band  strong  motion  instrument.  The  response  at  the  surface  for  finite  differences 
and  Cagniard  are  compared  in  Figure  10,  demonstrating  the  usefulness  of  finite 
differences  in  modeling  strong  motions. 

The  synthetics  for  small  ranges  have  not  been  included  in  either  Figure  9  or 
10,  because  the  line  to  point  source  mapping  approximation  breaks  down  at  near 
vertical  take-off  angles,  see  Appendix  A.  The  parameter  k,  used  in  the  mapping, 
is  fixed  at  0.6  throughout  this  investigation,  because  it  seems  to  give  good  agree¬ 
ment  to  other  synthetic  methods  (Figures  9  and  10)  at  the  take-off  angles  of 
interest  in  this  study. 

Canonical  basin  models 

Having  demonstrated,  briefly,  the  accuracy  of  the  finite-difference  and  two- 
dimensional  Kirchhoff  methods,  it  is  now  appropriate  to  examine  some  canonical 
models  of  basins,  to  understand  what  effects  various  geometries  may  have.  Four 
models  of  basin  boundaries  are  shown  in  Figure  11;  they  are  variations  of  a  layer 
over  a  half  space.  The  half  space  parameters  are  p-velocity  4.6  km/s,  s-velocity 
2.7  km/s,  and  density  2.7  g/cc;  and  are  representative  of  crystalline  rock.  The 
layer  parameters  are  chosen  to  give  a  p-velocity  ratio  of  one  to  two,  with  the 
other  parameters  appropriate  for  a  sediment  with  that  p-velocity;  p-velocity  2.3 
km/s,  s-velocity  1.1  km/s,  density  2.0  g/cc.  The  source  is  fixed  3.75  km  from  the 
basin  boundary  at  a  depth  of  0.75  km  in  each  case.  Surface  strong  motions  are 
computed  at  one  km  intervals  from  the  source  to  the  edge  of  the  model  and  com¬ 
pared  with  the  uniform  layer  over  a  half  space  model.  Examination  of  the  strong 
motions  in  Figure  12  shows  some  of  the  differences  between  the  four  basin  termi- 
nation  models.  Both  of  the  sharper  basin  terminations  pass  more  energy  across 
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the  termination;  the  gradual  terminations  allow  waves  of  only  about  half  the 
amplitude  to  pass  across.  Nevertheless,  all  the  models  of  basin  termination  cause 
large  drops  in  amplitude  as  the  wavefronts  cross  the  boundary.  Part  of  this 
energy  is  reflected  back  across  the  basin  (this  energy  cannot  be  corrected  to  three 
dimensions  properly  and  will  have  higher  than  its  true  amplitude,  especially  as  it 
approaches  the  source  position),  but  much  of  it  is  scattered  to  teleseismic  dis¬ 
tances.  The  surface  waves  are  not  well-developed  here  because  there  is  no  low- 
velocity  surface  layer  to  reduce  the  direct  wave  and  enhance  amplitude  and  dura¬ 
tion  of  the  Rayleigh  wave.  The  primary  conclusion  to  be  drawn  from  the  strong 
motions  is  that  gradual  basin  terminations  will  have  the  greatest  effect  on  energy 
crossing  that  boundary,  but  that  all  large  contrasts  across  great  changes  in  basin 
depth  will  cause  large  reductions  in  transmitted  strong  motion  amplitudes. 

Now,  using  the  two-dimensional  Kirchhoff  technique  derived  above,  the  scat¬ 
tered  energy  may  be  examined  at  teleseismic  distances.  Figure  13  shows  the 
teleseismic  p-wave  seismograms  computed  at  a  distance  through  the  half-space  of 
1000  km  from  the  source.  These  displacement  records  include  the  RDP  source,  a 
WWSSN  short-period  instrument  response  and  attenuation  with  T  *  of  one.  The 
seismograms  are  all  normalized  to  the  flat-layer  response  (ie:  the  response  if  the 
model  at  the  source  were  extended  laterally  without  termination).  The  distance 
of  1000  km  is  sufficient  for  these  half-space  teleseismic  calculations  because  it  is 
two  orders  of  magnitude  greater  than  the  length  of  the  integration  contour,  and 
therefore  the  waveform  will  no  longer  change  with  distance,  only  the  amplitude 
will  change  by  geometric  spreading.  This  will  be  discussed  more  fully  later. 
Differences  in  the  synthetics  for  the  various  models  consist  primarily  of  small 
amplitude  changes  at  the  frequencies  involved  here.  Even  broad  band  responses 
do  not  show  large  changes  in  teleseismic  records  (Figure  14).  This  indicates  that 
at  teleseismic  distances,  for  explosions,  velocity  contrast  at  the  boundary  and  the 
overall  dimensions  of  the  boundary  are  of  primary  importance,  not  the  precise 
shape  of  the  boundary.  On  the  other  hand,  increasing  the  amplitude  and  dura¬ 
tion  of  the  surface  wave  before  it  interacts  with  the  boundary  increases  the 
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overall  effect  on  the  teleseismic  waveform  and  increases  the  variation  with  boun¬ 
dary  geometry,  as  we  will  show  later. 

Some  of  this  scattered  energy  is  also  found  as  teleseismic  sv-waves.  Figure 
15  shows  the  s-wave  responses  at  the  same  teleseismic  distances.  These  waves  are 
of  the  same  order  of  magnitude  in  amplitude  as  the  scattered  p-waves.  Also, 
because  they  are  not  dominated  by  P  and  pP  (pS  is  small  at  near  vertical  take-off 
angles),  they  show  a  great  deal  more  variation  with  the  type  of  boundary  chosen. 
Such  high  amplitude  scattered  sv-waves  could  complicate  the  analysis  of  sources 
such  as  earthquakes  which  generate  direct  s-waves,  but  these  waves  could  be 
important  when  studying  explosions,  which  can  only  generate  sv-waves  through 
structural  interaction. 

In  demonstration  of  how  two-dimensional  Kirchhoff  constructs  the  telese- 
israic  waveform,  Figure  16  shows  a  series  of  synthetic  velocity  seismograms  at 
four  take-off  angles  (5,  10,  15  and  20  degrees)  for  five  cases  in  which  more  of  the 
line  elements  along  the  contour  are  included  in  each  subsequent  calculation. 
That  is,  case  0  includes  contributions  from  the  single  line  element  at  17  degrees 
take-off  angle.  Case  1  includes  the  two  nearest  neighbors,  case  2  includes  11  ele¬ 
ments,  case  3  includes  25,  and  case  4  includes  all  elements.  Up  to  case  2,  the  P- 
pP  arrivals  increase  by  addition  of  energy  from  a  narrow  range  of  angles  provid¬ 
ing  energy  propagating  away  from  the  very  near  vicinity  of  the  source.  Energy 
scattered  from  structure  is  not  propagating  in  the  direction  of  the  receivers  over 
this  portion  of  the  contour.  The  elements  used  in  case  2  should  include  most  the 
energy  for  the  P-pP  arrival,  and  remaining  contributions  to  that  phase  are  small, 
as  seen  by  comparing  peak  amplitudes  in  the  remaining  cases.  Case  3  is  the  first 
to  include  significant  energy  from  scattering.  Notice  that  the  contour  elements 
included  in  case  three  contain  all  contributions  from  the  basin  boundary  up  to  20 
degrees  take-off  angle.  For  this  reason,  the  contributions  from  the  remaining  ele¬ 
ments  are  small  and  the  records  for  case  4  differ  little  from  those  of  case  3. 
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Figure  17  demonstrates  the  effect  of  distance  on  the  resulting  seismograms. 
At  the  smaller  distances  of  5  and  10  km,  cases  0  and  1,  the  contribution  from 
scattering  propagates  at  a  greatly  different  angle  than  the  energy  from  the  source. 
For  receivers  below  the  source,  the  scattered  energy  is  small  because  it  is  pro¬ 
pagating  backward  off  the  scatterer.  Further  from  the  source,  at  distances  of  100 
and  1000  km,  cases  2  and  3,  there  is  little  difference  between  records  at  different 
distances  because  the  take-off  angles  from  the  source  and  the  scatterer  become 
indistinguishable.  Source  and  scatterer  become  one.  These  synthetic  seismo¬ 
grams  contain  only  the  p-waves;  if  the  s-waves  were  included,  the  records  at  5 
and  10  km  would  become  very  complicated  and  would  differ  greatly  from  those  at 
teleseismic  distances. 

Application  -  Events  at  Yucca  Flat,  NTS 

Since  near  source  structural  interaction  causes  variations  in  teleseismic 
waveforms  due  to  surface  wave  interaction,  records  of  this  scattering  are  likely  to 
be  quite  common.  Yucca  Flat  at  NTS,  Nevada  is  one  example  of  a  structure  that 
may  generate  these  effects.  The  lateral  variation  in  the  Yucca  Flat  area  is  due  to 
the  basin  structure  of  the  valley.  Sources  at  NTS  are  typically  located  near  the 
center  of  the  valley,  three  to  five  kilometers  from  the  basin  termination.  As 
shown  in  the  canonical  models  of  basins  above,  this  situation  will  cause  strong 
interaction  of  the  surface  wave  with  the  basin  boundary.  At  the  surface,  the 
amplitude  of  the  surface  wave  will  be  reduced  and  the  frequency  content  will 
change  as  the  wave  encounters  hard  rock. 

As  discussed  above,  the  two-dimensional  methods  used  in  this  investigation 
produce  a  three-dimensional  response  from  point  sources.  The  conversion  derived 
in  Appendix  A  requires  that  the  source  used  to  drive  the  finite-difference  scheme 
include  higher  order  (non-isotropic)  line  source  terms  [Vidale  and  Helmberger, 
1986].  For  an  explosion,  these  terms  have  radiation  patterns  such  that  the 
response  is  correct  only  from  half  the  source.  In  the  modeling  below,  the  struc¬ 
ture  in  the  invalid  direction  is  restricted  to  flat  layers.  The  result  of  this 
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modeling  is  a  structure  which  is  locally  cylindrical,  because  the  response  in  and 
out  of  the  plane  and  in  the  invalid  direction  is  flat  layered  and  has  no  effect  on 
the  seismograms  produced. 

Figure  18  is  a  map  of  Yucca  Flat  that  shows  the  outlines  of  the  hard  rock 
outcrops  which  approximately  mark  the  the  basin  boundary.  Also  displayed  are 
the  various  event  locations  and  a  network  of  strong  motion  recorders  for  one 
event  (FLASK).  Note  that  the  locally  cylindrical  geometry  is  appropriate  here, 
especially  in  the  direction  of  the  WWSSN  station  MAT,  Japan.  The  line  AB 
marks  the  location  of  a  cross-section  shown  in  Figure  19.  The  cross-section  is  a 
view  looking  south,  and  the  positions  of  various  sources  and  receivers  used  in  this 
study  are  shown.  The  dashed  line  is  the  integration  contour  for  the  two- 
dimensional  Kirchhoff.  There  are  three  materials  in  this  model:  low-velocity 
alluvium,  volcanic  tuffs,  and  hard  rock  (Cretaceous  granites  and  Paleozoic  rocks). 
The  large  differences  between  the  velocities  of  these  materials  is  what  makes 
lateral  variation  important. 

The  local  cylindrical  symmetry  of  the  resulting  model  is  shown  in  Figure  20. 
The  source  falls  on  the  axis  of  the  cylinder,  and  the  line  source  to  point  source 
conversion  is  valid  only  in  the  direction  of  the  basin  boundary.  The  seismograms 
generated  for  the  source  positions  indicated  are  responses  for  this  three- 
dimensional  geometry. 

Figure  3  shows  some  of  the  strong  motion  data  from  the  FLASK  event. 
Both  the  locations  of  the  receivers  and  the  location  of  the  event  are  shown  on  the 
map  of  Figure  18.  The  data  shown  in  the  figure  are  vertical  velocity  records  from 
stations  three  to  four  kilometers  from  the  event.  The  basin  becomes  shallow 
between  the  source  and  the  position  of  the  westernmost  stations,  creating  a 
geometry  where  the  stations  closest  to  the  source  see  a  distant  basin  boundary, 
while  those  to  the  west  see  a  basin  boundary  directly  below.  In  the  figure,  the 
differences  between  the  two  directions  is  readily  apparent.  To  the  west,  the 
amplitudes  of  surface  waves  are  reduced  by  roughly  half  with  respect  to  those  to 


the  east.  It  also  appears  that  the  longer  period  surface  waves  are  reduced  more 
than  shorter  period  waves.  It  will  be  shown  that  the  energy  lost  from  the  surface 
waves  is  converted  to  body  waves  and  is  found  on  teleseismic  records. 

Some  strong  motion  synthetics  for  a  source  at  position  1  are  shown  in  Figure 
21.  Notice  that  as  the  Rayleigh  wave  encounters  the  basin  boundary,  the  ampli¬ 
tude  becomes  roughly  half  and  the  longer  period  surface  waves  are  greatly 
reduced.  The  source  parameters  mentioned  in  the  figure  are  used  in  a  Helm- 
berger  and  Hadley  [1981]  reduced  displacement  potential.  This  RDP  source  is 
approximately  that  expected  for  FLASK  and  similar  events  in  Yucca  Flat,  and  is 
used  throughout  this  investigation  as  a  representative  source  description. 

Direct  comparisons  of  various  source-receiver  pairs  and  the  data  are  given  in 
Figure  22.  The  data  have  more  energy  in  the  higher  frequencies  due  to  limita¬ 
tions  of  the  finite-difference  method.  Finite  grid  spacing  results  in  a  maximum 
frequency  which  can  be  propagated  accurately.  Therefore,  the  source  used  is  con¬ 
volved  with  a  Gaussian  window  to  remove  higher  frequencies.  Nevertheless,  the 
absolute  amplitudes  and  many  of  the  primary  features  of  each  record  are  accu¬ 
rately  modeled.  The  record  for  station  795,  is  best  modeled  by  a  source  at  posi¬ 
tion  2  (out  of  positions  0  through  3  spaced  one  kilometer  apart)  at  a  distance  of 
three  kilometers.  This  is  one  of  the  better  waveform  fits  produced  by  the  model. 
Station  791  is  best  modeled  by  a  source  at  position  3  (closest  to  the  basin  boun¬ 
dary)  at  a  distance  of  3.5  kilometers.  Stations  791,  793  and  795  are  north  and 
east  of  the  event.  This  is  the  range  of  directions,  as  argued  from  the  map  of  Fig¬ 
ure  18,  which  is  most  like  a  locally  cylindrical  geometry.  Station  789  is  to  the 
south.  Here,  the  basin  boundary  is  further  from  the  source.  Thus,  source  posi¬ 
tion  0  gives  the  most  accurate  result. 

Now  the  two-dimensional  Kirchhoff  may  be  applied  to  see  if  this  model  is 
consistent  also  with  the  teleseismic  data.  Figure  2  shows  the  data  set  to  be  con¬ 
sidered  here.  These  are  short  period  WWSSN  records  recorded  at  station  MAT, 
Japan.  Seven  different  events  from  various  sites  within  the  valley  are  shown. 


Notice  the  additional  energy  arriving  after  P-pP.  This  energy  is  present  at  very 
large  amplitudes  on  all  records  except  PORTMANTEAU.  PORTMANTEAU 
occurred  at  the  basin  boundary  to  the  northeast  and  therefore  may  be  considered 
primarily  a  flat  layer  case.  Note  also  that  while  P-pP  is  similar  for  all  the  events, 
the  later  phases  vary  substantially  as  a  function  of  position.  Yet  when  sources 
are  close  together  (KEELSON  and  OSCURO),  differences  are  small. 

Figure  23  shows  some  of  the  two-dimensional  Kirchhoff  results.  These  are 
point  source  displacement  records  and  have  been  convolved  with  a  WWSSN  short 
period  instrument  response  and  a  Q-operator  with  T*=l.  Take  off  angles  of  15° 
and  20°  were  chosen  to  bracket  that  appropriate  for  MAT  (17.8°).  All  amplitudes 
are  relative  to  the  flat  layer  case  at  a  take-off  angle  of  0°  (the  flat  layer  record 
shown  is  at  15°).  The  important  observation  here  is  that  the  records  vary  far 
more  by  moving  the  source  one  kilometer  within  the  basin  than  by  changing  the 
take-off  angle  by  five  degrees.  The  energy  which  causes  the  waveform  variations 
comes  from  the  conversion  of  surface  wave  energy  at  the  basin  boundary.  Figure 
24  shows  a  direct  comparison  of  the  data  and  model.  KEELSON  and  OSCURO 
are  modeled  accurately  by  the  same  synthetic  record,  demonstrating  repeatabil¬ 
ity,  and  the  flat  layer  comparison  with  PORTMANTEAU  is  good  as  expected. 

Discussion  and  Conclusions 

This  paper  presents  a  new  method  for  calculating  teleseismic  waveforms, 
using  an  interface  of  numerical  and  analytic  computations.  As  presented,  the 
two-dimensional  Kirchhoff  method  relies  entirely  on  local  structural  effects  to 
change  the  teleseismic  waveform.  This  method  can  be  made  computationally  fast 
and  efficient.  Using  anisotropic  sources  and  filtering  the  resulting  seismograms 
produce  three  dimensional,  point  source  records.  The  finite-difference  calculation 
produces  both  surface  strong  motion  records  and  the  functions  needed  for  the 
teleseismic  synthetics.  In  this  way,  a  model  of  local  structure  can  be  compared 
with  both  local  and  teleseismic  data  simultaneously. 


In  this  paper,  we  concentrated  on  the  local  structural  effects  in  modeling 
teleseismic  waveforms.  More  realistic  path  effects  can  be  included  in  the  calcula¬ 
tions  by  modifying  Equation  39  using  a  convolution  with  an  appropriate  kernel 
computed  by  various  analytical  ray  techniques.  Such  interfacing  would  be  useful, 
for  example,  in  modeling  the  complexities  observed  by  Lynnes  and  Lay  [1987] 
(also  on  this  issue).  This  problem  will  be  addressed  in  future  publications. 

The  case  for  using  only  local  structure  in  the  teleseismic  calculation  is  strong 
for  Yucca  Flat  because  the  local  strong  motion  and  teleseismic  data  are  modeled 
simultaneously.  Only  the  local  structure  need  be  considered  to  make  a  good  syn¬ 
thetic  strong  motion  seismogram.  If  a  local  structure  models  well  both  the  strong 
motion  records  and  the  teleseismic  records  then  the  effect  of  local  structure  on 
the  teleseismic  waveform  is  substantiated.  A  relatively  simple  model  of  the  basin 
at  Yucca  Flat  produces  synthetic  seismograms  which  match  very  well  those 
observed  at  teleseismic  and  local  distances.  Several  features  of  these  records  are 
explained.  The  variation  of  surface  wave  amplitude  with  azimuth  is  shown  to  be 
the  result  of  interaction  with  the  basin  boundary.  A  drop  in  amplitude  across 
the  basin  boundary  occurs  as  surface  wave  energy  is  converted  to  body  waves. 
These  scattered  phases  are  observed  on  teleseismic  records  shortly  after  P-pP. 
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Appendix  A:  Analytic  source  function  for  2d  finite  difference 


Asymptotic  source  theory  produces  p-wave  point  source  displacement  potentials 
in  three  dimensions  of  the  following  forms,  assuming  a  step  function  source 
[Helmberger,  1983] 


(A.1) 


and  for  a  delta  function  source 
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In  these  equations,  p  is  given  by 


f  -  £  + 


(A.3) 


and  1 1  is  given  by 


(A.4) 


The  conventions  are  x  positive  to  the  right,  «  positive  to  the  right,  *  positive 
down,  but  w  positive  up.  The  variables  are:  x,  u  horizontal  distance  and  dis¬ 
placement;  i,  w  vertical  distance  and  displacement;  a  p-wave  velocity;  4>  poten¬ 
tial;  and  R  =(x2+*2)l/a. 


The  method  used  in  this  paper  requires  the  delta  function  response  in  dis¬ 
placement,  which  is  obtained  by  taking  derivatives  of  the  potential. 


(A.5) 


To  simplify  the  taking  of  these  derivatives,  first  take  the  Laplace  transform  of 
Equation  A.2. 
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Now,  taking  the  derivatives  A.5  of  Equation  A.6  and  keeping  only  the  far-field 


terms 
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(A.7) 


where  <  is  the  signum  function. 
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Returning  A.7  to  the  time  domain, 
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(A.8) 


Part  of  Equation  A.2  can  be  evaluated  before  assembling  the  displacements,  that 


t\l± 

l  *j  dt 


H(t  -R  /a) 
/c)5 


(A.0) 


A  final  stipulation  is  that  one  derivative  with  respect  to  time  will  be  applied  later 
to  the  source  time  function  before  convolving  it  with  the  result,  so  one  derivative 
can  be  dropped  here.  Then  using  A.2  and  A.9  in  Equation  A.8: 
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Now  let 
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Then,  combining  A. 10  and  A.ll, 
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Approximate  x/jT  as 
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where  k  is  a  parameter  that  allows  a  better  fit  to  sfp  over  selected  ranges  of 
take-off  angle.  Using  this  approximation,  set 


Re  (p  sfp  )  =  Re  p 
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(A.15) 


Similarly, 
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If  I  is  the  take-off  angle  (the  direction  of  Ft),  measured  counter-clockwise  from 
positive  z,  then  p0— *in(i)/a.  Take  p0=  l/o  for  « — >r/2  in  Equations  A.15  and  A.16, 
and  apply  these  to  Equation  A.  13, 
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Let  (R  /at  )2—Ta  and  |z|<==*,then 
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It  is  instructive  at  this  point  to  do  dimensional  analysis  on  Equation  A.17. 
Let  units  of  length  be  denoted  /  and  units  of  time  t .  Then, 
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Convert  this  result  to  point  source  form  using  Equation  A.12. 

(  111  s/T  A 

'*>’*  “  7T  t7T  s/P 

=  -jr  (A.  19) 

Now  convolve  the  result  with  the  correct  time  function,  the  reduced  displacement 
potential  (RDP)  function.  The  definition  of  the  RDP  is 


<KR  -0  =  where  0(f)  =  0«(l-e '“)(l+*f  +...) 

where  0,*  has  dimensions  of  volume.  From  Equation  A.5,  displacement  is 

94  0  1  dr{> 

*  ~  9R  ~  R*  +  R  a  dt 

Applying  this  to  the  derivation  above,  true  displacement  is  expressed  as 


«(/?,<)  =  -V<0*- 


(A.20) 


(A.21) 


(A.22) 


where  the  derivative  is  equivalent  to  the  expression  of  Equation  A.19.  Continu¬ 
ing  with  the  dimensional  analysis,  Equations  A.19  and  A.22  give  the  result 


/s4-f  -  U 


(A.23) 


But,  in  the  derivation  of  Equations  A. 10  and  A. 11,  one  time  derivative  was 
dropped  in  order  that  it  could  later  be  applied  to  the  source  time  function. 
Replacing  4  with  i  4/dt  in  equation  A.22,  Equation  A.23  becomes 


t  l 2 


(A. 24) 


which  is  correct  for  displacement.  Instead  of  using  the  RDP,  a  moment  could  be 
applied  to  the  filtered  result.  Replace  0(f)  with  A/0l(r10/4Waa,  where  M0  is 
moment  in  dyne-cm,  d  is  density  in  g/cc,  and  a  is  p-velocity  in  km/s  for  dis¬ 
placement  in  cm. 


When  applying  this  source  to  the  finite  difference  method,  the  full  form  of 
the  source  is  assembled  in  stages.  The  form  of  the  source  entered  into  the  finite 
difference  grid  is  given  in  Equation  A.17.  The  value  of  k  is  chosen  at  the  start  of 
the  run  to  best  fit  the  required  range  of  take-off  angles.  Of  course,  since  the 
parameter  k  governs  a  linear  combination  between  two  source  terms,  a  line  iso¬ 
tropic  explosion  and  a  line  dipole  force,  the  results  of  two  runs  differing  only  in  k 
may  be  combined  to  give  any  desired  k.  Since  Equation  A.17  is  singular  at 
Ta—l,  that  is  t  =R  /o,  the  corresponding  point  for  the  numerical  time  series  must 
be  found  by  integrating  the  equation  and  matching  area.  This  singular  nature 
also  introduces  energy  at  frequencies  too  high  for  the  grid  to  propagate.  To 
correct  this,  the  numerical  time  series  is  convolved  with  a  Gaussian  filter  before  it 
is  propagated. 

The  finite  difference  results  must  be  subjected  to  a  line  source  to  point 
source  conversion  filter.  This  filter  is  given  by  Equation  A.12.  It  is  restated  here 
where  C/liM  and  Ufaial  are  line  source  and  point  source  wavefields. 
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This  may  not  be  the  best  teleseismic  form,  considering  the  approximation  (Equa¬ 
tion  A.14)  of  v7.  Starting  with  Equation  A.2,  use  Equation  A.9  and  rearrange 
terms  to  obtain 


<j>j>  =  ^  Re  ( \fp  ) 
As  t  approaches  R  /a, 
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Use  Equation  A.14  with  k  =0.5  and  the  approximations 


to  get  the  form 
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2  7xR 


H(t-R  /a) 


The  reference  take-off  angle  i0  is  along  positive  * ,  that  is,  sin  =  1.  In  this  case, 
the  final  form  is 


4>p  —  ‘  )H (*  -R  / a)  (A.26) 

An  ideal  line  to  point  filter  would  cause  in  Equation  A.26  to  be  H(t-R  /o)/R 
for  all  angles  «  between  0°  and  180°.  A  better  filter  can  be  found  by  inspection. 
Replace  \fz  in  Equation  A.26  as  follows. 

4>?  —  2^/JT  +  *'n  ,  )^(<~^/0t)  (A.27) 

Then  for  both  i'-=0°  and  i— =90°  it  is  immediately  found  that  is  H(t-R  /a)/R . 
Finally,  applying  the  substitution  used  to  obtain  Equation  A.27  to  the  filter  given 
in  Equation  A.25, 
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Figure  1.  Schematic  diagram  displaying  energy  paths  lor  a)  flat  layered 
model  versus  b)  laterally  varying  structure.  The  model  is  two-dimensional. 
The  figure  demonstrates  the  motivation  for  developing  a  two-dimensional 
teleseismic  method  that  accounts  for  local  structural  variations. 

Figure  2.  Teleseismic  data  at  MAT.  These  are  WWSSN  short  period  vert¬ 
ical  seismograms  recorded  at  station  MAT,  Japan  for  seven  different 
events  at  various  locations  within  Yucca  Flat  (see  Figure  18).  There  are 
very  strong  secondary  phases  appearing  after  P-pP.  While  P-pP  seems 
relatively  consistent  among  the  records,  the  secondary  phases  are  not,  indi¬ 
cating  the  presence  of  scattering  structure  near  the  source. 

Figure  3.  Strong  motion  data  for  FLASK.  Only  the  vertical  velocity 
records  are  shown  here.  The  stations  (numbered  stations  ranging  from  781 
to  795)  are  each  marked  as  an  x  on  the  map  of  Figure  18.  The  source  is  in 
the  center.  Such  data  strongly  argue  for  lateral  variation  because  of  the 
strong  lateral  contrast  in  peak  amplitude  and  duration  of  the  Rayleigh 
wave. 

Figure  4.  Green’s  transformation  for  P  not  in  domain  D.  The  contour  T  is 
chosen  such  that  there  are  no  singularities  of  the  integrand  in  the  area  D 
that  T  encloses.  The  variable  r  in  the  transformation  is  the  distance  from 
P  to  any  point  on  T  or  in  D.  The  vector  n  is  the  outward  directed  normal 

to  r. 

Figure  5.  Green’s  transformation  for  P  in  D.  Here,  the  contour  a  is  intro¬ 
duced  in  addition  to  T  around  the  point  P.  Now  the  domain  D’  is  that 
area  enclosed  between  contour  T  and  contour  a.  These  contours  are  freely 
deformable  and  a  is  chosen  to  be  circular  with  radius  t  centered  on  point 
P.  The  limit  as  t  approaches  0  will  produce  the  desired  integral  formula. 
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Figure  6.  Geometry  for  three-dimensional  analytic  evaluation  of  the  Kir- 
chhoff  integral.  The  interface  surface  is  S,  with  solid  angle  ft  tracing  out 
dS.  The  normal  to  the  surface  is  n.  The  source,  at  distance  r0  from  S, 
generates  potential  4>  which  produces  signal  <£(S)  on  the  surface.  Kirchhoff 
integration  gives  the  result  $(P)  at  point  P,  at  a  distance  r  from  S. 

Figure  7.  Geometry  for  two-dimensional  analytic  evaluation  of  the  two- 
dimensional  Kirchhoff  integral.  The  line  interface  is  T,  with  angle  6  trac¬ 
ing  out  dl.  The  normal  to  the  contour  is  n.  The  source,  at  distance  r0 
from  T,  generates  potential  4>  which  produces  signal  4>{T)  on  the  surface. 
Two-dimensional  Kirchhoff  integration  gives  the  result  $(P)  at  point  P,  at 
a  distance  r  from  T. 

Figure  8.  The  geometry  for  two-dimensional  Kirchhoff  along  a  flat  inter¬ 
face  in  a  half  space.  The  contour  F  has  been  deformed  such  that  it  follows 
a  line  rt  parallel  to  the  x  axis  from  x=-oo  to  x=+oo  and  is  closed  at 
P  =oo  (r2)  in  the  direction  of  z=+oo.  The  integrand  becomes  trivial 
along  r2,  reducing  the  integral  to  an  infinite  definite  integral  along  x.  This 
is  further  reduced  for  numerical  application  to  a  finite  sum  along  x. 

Figure  9.  Test  results  for  Lamb’s  problem.  Three  methods  are  compared 
in  this  figure:  two-dimensional  Kirchhoff,  finite  differences,  and  Cagniard  - 
de  Hoop.  The  records  are  responses  in  a  half  space  with  a  source  1  km 
deep  for  deeply  buried  receivers  (6  km)  at  the  range  of  horizontal  distances 
indicated.  The  velocities  of  the  medium  are  4.0  km/s  p-wave  and  2.3 
km/s  s-wave  with  a  density  of  2.7  gm/cc.  The  source  is  an  explosion  with 
RDP  parameters  K=12.0,  B=1.0  and  0^=10'°. 

Figure  10.  Test  results  for  Lamb’s  problem.  The  medium  is  the  same  as 

that  of  Figure  9.  The  receivers  are  on  the  surface  and  only  finite 

% 

differences  and  Cagniard  are  compared.  This  demonstrates  the  accuracy 
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of  finite-difference  strong  motions. 


Figure  11.  Four  models  of  basin  boundaries.  These  are  the  four  canonical 
models  which  are  used  to  demonstrate  the  effect  of  various  basin  termina¬ 
tions.  The  structure,  apart  from  the  boundaries,  is  a  layer  over  a  half 
space.  The  layer  has  a  p-velocity  2.3  km/s,  s- velocity  1.1  km/s  and  den¬ 
sity  2.0  gm/cc.  The  half  space  has  a  p-velocity  4.6  km/s,  s- velocity  2.7 
km/s  and  density  2.7  gm/cc.  The  star  represents  the  position  of  the 
source,  always  at  a  depth  of  750  m.  The  inverted  triangles  are  the  posi¬ 
tions  of  strong  motion  instruments.  The  broken  line  is  the  contour  along 
which  two-dimensional  Kirchhoff  is  performed.  The  distances  labeled  on 
model  1  are  in  km. 

Figure  12a.  Strong  motion  results  for  canonical  basin  models.  These 
seismograms  are  velocity  records  at  the  surface  at  distances  1,  3,  5  and  7 
km  from  the  source,  for  the  first  two  models  shown  in  Figure  11.  An  RDP 
source  with  K— 12,  B=l,  and  ^'oo=1010  has  been  used.  The  number  to 
the  right  of  each  trace  is  the  maximum  amplitude  in  cm/s  along  that 
trace.  Each  trace  from  each  model  (heavy  line)  is  compared  directly  to  the 
corresponding  trace  for  flat  layers  (light  line).  Two  important  observations 
are  made.  First,  the  amplitude  of  the  surface  wave  drops  abruptly  at  the 
basin  boundary,  and  second,  this  drop  is  greatest  for  the  least  dipping 
boundary. 

Figure  12b.  Strong  motion  results  for  basin  models.  This  figure  is  the 
same  as  Figure  12a,  but  compares  the  remaining  two  models  shown  in  Fig¬ 
ure  11. 

Figure  13.  Teleseismic  p-wave  results  for  canonical  basin  models.  These 
seismograms  are  displacement  records  convolved  with  a  WWSSN  short 
period  instrument  response  and  include  the  same  RDP  source  used  in 
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Figure  12.  They  are  also  convolved  with  a  Q  operator  with  T*=l.  They 
are  at  a  ray  length  of  1000  km  from  the  source,  and  at  take-off  angles  10°, 
15°  and  20°.  The  number  to  the  right  of  each  record  is  the  peak  ampli¬ 
tude  normalized  to  the  peak  amplitude  of  a  flat  layer  model  at  0°  take  of 
angle.  The  simple  appearance  of  these  seismograms  is  misleading,  they 
vary  significantly  from  the  flat  layer  response  shown  as  the  fifth  record  in 
each  series.  The  first  four  records  are  from  each  model,  in  order,  as  shown 
in  Figure  11. 

Figure  14.  Teleseismic  p-wave  results  for  basin  models.  This  is  identical  to 
Figure  13,  except  that  no  Q  or  instrument  response  has  been  convolved 
with  the  record.  This  makes  the  records  broad  band.  Differences  between 
the  records,  especially  secondary  arrivals,  are  now  more  obvious. 

Figure  15.  Teleseismic  s-wave  results  for  basin  models.  The  purpose  of 
this  figure  is  to  demonstrate  the  production  of  teleseismic  s-waves, 
although  the  true  character  of  teleseismic  s-waves  observed  is  not  accu¬ 
rately  reflected  here.  The  figure  is  set  up  the  same  as  Figure  13,  but  the 
value  of  T*  is  far  too  low  for  the  earth  for  s-waves.  The  sv  arrivals  are 
partly  pS,  but  this  phase  is  small  at  small  take-off  angles.  Much  of  the 
energy  shown  is  the  result  of  basin  boundary  interactions. 

Figure  16.  Contributions  of  various  line  elements  to  the  final  seismogram. 
The  model  is  model  4  from  Figure  11.  Case  0  is  the  result  of  the  contribu¬ 
tion  from  one  line  element  at  17°.  Each  case  adds  more  line  elements  to 
the  final  seismogram.  Cases  0  through  2  progressively  build  the  P-pP 
arrival;  case  3  adds  the  energy  from  scattering.  The  records  in  all  cases  are 
broad  band  velocity  records  at  5,  10,  15  and  20  degrees  take-off  angles. 
They  include  an  RDP  source  as  described  with  Figure  12. 
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Figure  17.  Effect  of  distance.  Cases  0  through  3  represent  progressively 
greater  distance  from  the  source  (5,  10,  100  and  1000  km).  These  are 
broad  band  displacement  records  at  5,  10,  15  and  20  degrees  take-off  angle. 
Cases  3  and  4  are  similar  because  at  great  distance,  thus  the  source 
becomes  indistinguishable  from  the  scatterer  in  terms  of  take-off  angle. 


Figure  18.  Map  of  Yucca  Flat.  The  contours  show  the  borders  of  hard 
rock  outcrops.  The  filled  circles  are  the  locations  of  several  events  listed 
adjacent  to  the  figure.  Each  strong  motion  station  used  to  record  the 
FLASK  event  (number  12)  is  represented  by  an  x.  The  arc  indicates  the 
general  shape  of  the  basin  boundary  as  seen  from  above.  The  arrow  indi¬ 


cates  the  direction  of  WWSSN  station  MAT,  Japan. 


Figure  19.  Cross  section  of  Yucca  Flat.  The  cross  section  shown  here  is  a 


generalization  of  the  known  structure  at  Yucca  Flat.  The  three  layers 
from  the  surface  downward  are  alluvium,  volcanic  tuff  and  hard  rock. 


Hard  rock  is  a  broad  term  used  here  to  describe  Mesozoic  granitic  rock  and 
Paleozoic  rocks,  all  of  which  have  similar  elastic  properties  (velocities  and 
density).  The  stars  indicate  four  different  source  positions  used  (numbered 
0  to  3),  all  at  a  depth  of  875  m.  The  inverted  triangles  indicate  the  posi¬ 
tion  of  the  strong  motion  records  generated  by  the  finite-difference  calcula¬ 


tions. 


Figure  20.  Local  cylindrical  symmetry.  The  structure  of  Yucca  Flat  is 
shown  here  as  the  locally  cylindrical  result  required  by  the  filter  used  to 
convert  line  to  point  responses.  Comparing  this  to  the  map  of  Figure  18, 
locally  cylindrical  structure  is  desirable  for  the  case  of  Yucca  Flat. 


Figure  21.  Synthetic  strong  motion  records  for  FLASK.  These  records, 
generated  by  the  finite- difference  method,  include  an  RDP  source  with 
K=12,  B=1  and  \poo=l010.  These  are  only  for  source  position  1  (of  0  to  3 
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on  Figure  19),  and  for  the  stations  shown.  It  is  important  to  see  that  the 
peak  amplitude  drops  sharply  across  the  boundary,  that  the  duration  of 
the  Rayleigh  wave  is  reduced,  and  that  the  Rayleigh  wave  appears  to  lose 
relatively  more  of  the  lower  frequencies  as  it  crosses  the  basin  boundary. 

Figure  22.  Direct  comparison  of  strong  ground  motion.  This  figure  shows 
some  direct  observed  to  synthetic  comparisons  of  strong  ground  motions 
for  the  FLASK  event.  This  demonstrates  that,  despite  the  simplicity  of 
the  model  used,  the  resulting  strong  motions  are  an  accurate  representa¬ 
tion  of  those  observed. 

Figure  23.  Synthetic  teleseismic  records  for  Yucca  Flat  events.  A  short 
period  WWSSN  response  and  a  Q  operator  with  T*=l  have  been  con¬ 
volved  into  the  records.  The  RDP  source  mentioned  in  Figure  21  is  also 
included.  The  peak  amplitudes  are  normalized  to  that  for  a  flat  layer 
response  at  a  take-off  angle  of  0°.  The  flat  layer  record  shown  is  at  15° 
The  important  observation  here  is  that  moving  the  source  1  km  within  the 
basin  is  far  more  important  than  changing  the  take-off  angle  5°. 

Figure  24.  Direct  comparison  of  teleseismic  records.  These  comparisons 
demonstrate  the  accuracy  of  the  new  method.  Although  the  p-waves  have 
passed  through  the  mantle  and  receiver  structure,  the  majority  of  the 
energy  in  the  observed  records  is  explained  by  the  near  source  structure  at 
Yucca  Flat.  The  repeatability  of  the  method  is  demonstrated  by  KEEL¬ 
SON  and  OSCURO.  These  events  are  located  close  together  and  are 
modeled  well  by  the  same  synthetic  record. 
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Abstract 


The  methods  for  Representation  Theorem  (RT)  coupling  of  finite  element  (FE) 
or  finite  difference  calculations  and  Harkrider’s  (Harkrider  1964,  1970)  propagator 
matrix  method  calculations  to  produce  a  hybrid  method  for  propagation  of  SH  mode 
sum  seismograms  across  paths  that  contain  regions  of  non  plane-layered  structure  are 
explained  and  developed.  The  coupling  methods  explained  in  detail  use  a  2-D  Carte¬ 
sian  FE  formulation.  Analogous  methods  for  the  3-D  method  follow  directly.  Exten¬ 
sive  tests  illustrating  the  validity  and  accuracy  of  the  implementation  of  these  cou¬ 
pling  methods  are  discussed.  These  hybrid  techniques  are  developed  to  study  the 
propagation  of  surface  waves  across  regional  transition  zones  or  other  heterogeneities 
that  exist  in  part  of  a  longer,  mostly  plane-layered,  path.  The  effects  of  a  thinning  or 
thickening  of  the  crustal  layer  on  the  propagation  of  Lg  mode  sum  seismograms  have 
been  examined  in  this  study.  The  thinning  or  thickening  of  the  crustal  layer  is  used 
as  a  simple  model  of  ocean  continent  transitions.  The  Lg  phase  is  of  particular 
interest  since  it  is  used  in  several  important  applications  such  as  mapping  the  extent 
of  continental  crust,  magnitude  determination,  and  discrimination  between  explosive 
and  earthquake  sources.  The  understanding  of  the  observations  that  Lg  wave  is 
attenuated  completely  when  the  propagation  path  includes  an  oceanic  portion  of 
length  greater  than  one  hundred  to  two  hundred  kilometers  or  a  region  of  complex 
crustal  structure  is  not  complete,  and  a  clear  explanation  of  these  phenomena  could 
have  important  consequences  for  all  these  types  of  studies.  The  transition  model  cal¬ 
culations  done  in  this  study  show  that  passage  through  a  region  of  thinning  crustal 
thickness,  the  model  for  a  continent  to  ocean  transition,  increases  the  amplitude  and 
coda  length  of  the  Lg  wave  at  the  surface,  and  allows  much  of  the  modal  energy 
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trapped  is  the  crust,  which  forms  the  L(  phase,  to  escape  into  the  subcrustal  layers 
as  body  waves  or  other  downgoing  phases.  The  magnitude  of  both  these  effects 
increases  as  the  length  of  the  transition  increases  or  the  slope  of  the  layer  boundaries 
decrease.  The  passage  of  the  wavefront  exiting  the  continent  to  ocean  transition 
region  through  the  oceanic  structure  allows  further  energy  to  escape  from  the  crustal 
layer,  and  produces  a  decrease  in  L(  amplitude  at  the  surface  as  the  length  of  the  oce¬ 
anic  path  increases.  The  amplitude  decrease  is  maximum  near  the  transition  region 
and  decreases  with  distance  from  it.  Passage  through  a  region  of  thickening  crust, 
the  model  of  a  ocean  to  continent  transition,  causes  a  rapid  decrease  in  the  Lf  ampli¬ 
tude  at  the  surface  of  the  crust.  The  energy  previously  trapped  in  the  oceanic  crustal 
layer  spreads  throughout  the  thickening  crustal  layer,  and  any  amplitude  which  has 
been  traveling  through  the  subcrustal  layer  but  has  not  reached  depths  below  the 
base  of  the  continental  crust  is  transmitted  back  into  the  continental  crust.  The 
attenuation  of  Lf  at  the  crustal  surface  along  a  partially  oceanic  path  occurs  in  the 
oceanic  structure  and  in  the  ocean  to  continent  transition  region.  The  attenuation 
at  the  surface  depends  in  part  on  the  escape  of  energy  at  depth  through  the  continent 
to  ocean  transition  region  into  the  underlaying  half-space.  The  total  attenuation  of 
Lg  due  to  propagation  through  a  forward  transition  followed  by  a  reverse  transition 
is  at  most  a  factor  of  four  to  six.  This  is  inadequate  to  explain  the  observed  attenua¬ 
tion  of  Lg.  Thus,  additional  effects,  other  than  geometry  must  be  considered  to  pro¬ 
vide  a  complete  explanation  of  the  attenuation  of  Lg. 
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Chapter  3 

Finite  Element  to  Modal  Propagator  Matrix  Coupling:  Testa  of  Accuracy  and 
Applications  to  the  Transmission  of  Lg  along  Partially  Oceanic  Paths 


Introduction 

In  this  chapter  the  effects  of  the  length  of  the  intermediate  path  between  the  con¬ 
tinent  ocean  and  ocean  continent  transition  regions  is  investigated.  First  the  results  of 
two  FE  calculations  with  different  intermediate  path  lengths  are  presented  and  com¬ 
pared.  These  examples  are  contrasted  with  the  path  length  used  in  the  previous 
chapter.  Then  the  RT  integration  method  is  discussed  and  explained.  First  analytic 
expressions  for  the  stress  components  of  double  couple  and  line  sources  are  derived, 
then  the  expressions  for  displacement  and  stress  line  source  Green’s  functions  are 
determined.  These  expressions  are  used  to  illustrate  the  validity  and  determine  the 
accuracy  of  the  RT  coupling  method.  The  RT  coupling  method  can  be  used  to  con¬ 
tinue  the  propagation  of  FE  results  through  a  layered  structure  using  the  displacement 
and  stress  Green’s  functions  for  the  remaining  path  length  and  the  displacement  and 
stress  time  histories  recorded  at  a  column  of  element  centers  within  the  FE  grid. 

FE  Results:  Effect  of  oceanic  path  length  on  the  attenuation  of  Lg 

Results  of  two  FE  calculations  including  both  forward  and  reverse  transition 
regions  and  an  intermediate  region  of  oceanic  structure  are  presented  here.  Many  of 
the  details  of  transmission  of  wavefield  through  a  forward  transition  region  or  a 
reverse  transition  region  have  previously  been  discussed  and  will  not  be  repeated  here. 
The  vertical  extent  of  the  grid  and  thus  that  of  each  time  slice  illustrated  for  these 
calculations  is  larger  than  that  for  the  results  discussed  in  the  previous  chapter.  This 
increase  in  vertical  extent  makes  the  disturbances  moving  down  through  the 
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underlying  half-space  easier  to  see.  The  change  in  the  angle  at  which  these  disturb¬ 
ances  travel  down  toward  the  base  of  the  grid  as  the  width  of  the  triangular  regions  of 
maxima  widen  later  in  the  wavefield  is  also  more  clearly  visible  in  these  examples. 

In  Figures  1  to  4  time  slices  recorded  during  a  calculation  using  a  grid  which 
included  fifty  kilometer  long  forward  and  reverse  transition  regions  separated  by  an 
intermediate  path,  whose  length  is  thirty  two  kilometers,  of  oceanic  structure.  In  Fig¬ 
ures  5  to  7  the  results  of  a  calculation  using  a  grid  containing  the  same  fifty  kilometer 
forward  and  reverse  transition  regions  separated  by  an  intermediate  path  of  seventy 
nine  kilometers  are  illustrated.  Figure  8  shows  the  RMS  amplitude  as  a  function  of 
distance  along  the  surface  of  each  of  these  models.  The  upper  plot  shows  the  ampli¬ 
tudes  for  the  thirty  two  kilometer  oceanic  path  and  the  lower  plot  shows  the  ampli¬ 
tudes  for  the  seventy  nine  kilometer  long  oceanic  path.  The  amplitude,  or  vertical 
scale  on  each  of  these  plots  is  identical,  however,  the  horizontal  scale,  measuring  dis¬ 
tance  along  the  surface  of  the  crustal  layer,  is  not  uniform.  The  vertical  bars  on  the 
plots  indicate  the  ends  of  the  regions  of  oceanic  structure.  The  horizontal  scales  are 
uniform  between  these  bars  and  uniform  but  different  outside  these  bars.  Between 
these  bars  the  horizontal  distance  scale  is  defined  so  that  the  total  length  of  the  oce¬ 
anic  path  in  the  illustrated  model  scales  exactly  to  the  distance  between  the  bars.  The 
scale  outside  the  area  defined  by  these  bars  is  such  that  the  distance  between  each 
point  in  the  regions  of  equally  spaced  points  is  ten  nodes. 

Examining  the  time  slices  for  the  two  examples  with  different  oceanic  path 
lengths  shows  that  the  major  difference  between  them  appears  to  be  the  amount  of 
energy  which  is  escaping  out  through  the  bottom  of  the  grid.  In  the  case  of  the 
shorter  path  length  only  the  disturbances  traveling  through  the  half-space  with  the 
steepest  angles  are  able  to  pass  out  of  the  system.  All  other  disturbances  are  largely 
or  completely  transmitted  back  into  the  continental  crustal  layer  since  they  have  not 
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propagated  downwards  far  enough  to  avoid  being  incident  of  the  mantle  crust  boun¬ 
dary  in  the  reverse  transition.  This  observation  provides  an  explanation  of  why  there 
is  a  critical  length  of  intermediate  oceanic  path  beyond  which  the  Lg  wave  does  not 
travel.  As  the  path  length  in  the  oceanic  region  increases,  disturbances  traveling 
through  the  crust  mantle  boundary  in  the  forward  transition  region  at  progressively 
shallower  angles  are  able  to  propagate  downwards  beyond  the  base  of  the  continental 
crustal  layer  and  escape  into  the  mantle.  This  means  that  less  energy  is  available  to 
be  transmitted  back  through  the  mantle  crust  boundary  of  the  reverse  transition  and 
reconverted  to  Lg  wave  energy.  This  trend  is  seen,  although  less  clearly,  in  the  for¬ 
ward  and  reverse  time  slices  in  chapter  2.  The  intermediate  oceanic  path  length  for 
that  example  is  one  hundred  sixty  nine  kilometers.  Disturbances  that  are  seen  to 
propagate  out  of  the  grid  in  that  long  path  example  can  be  traced,  even  in  the  exam¬ 
ple  with  the  longer  path  length  presented  here,  across  the  oceanic  path  length  back 
into  the  crust  through  the  reverse  transition  region.  This  implies  that  as  the  path 
length  increases  more  and  more  converted  energy  is  able  to  leave  the  system. 

There  is  another  contribution  of  the  oceanic  path  length  to  the  attenuation  of  Lg 
passing  through  an  oceanic  structure.  The  energy  that  is  transmitted  into  the  oceanic 
crustal  layer  is  not  all  in  the  form  of  modes  compatible  with  trapped  modes  in  the  oce¬ 
anic  crust.  Therefore,  as  the  wavefront  propagates  through  the  oceanic  structure 
further  energy  leaks  from  the  crustal  layer  and  propagates  downward  to  eventually 
escape  from  the  system. 

Calculation  of  Stress  Components 

To  test  the  RT  algorithm  for  coupling  FE  results  to  propagator  calculations  it  is 
first  necessary  to  determine  expressions  for  the  stress  seismograms  and  the  stress  and 
displacement  Green’s  functions  used  in  the  convolution  integral.  These  expressions  are 
developed  and  demonstrated  in  this  section.  The  use  of  these  expressions  in  the  RT 
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Figure  8:RMS  surface  amplitudes  for  the  two  FE  models  including  both  forward  and  reverse 
transitions.  The  upper  plot  shows  the  RMS  surface  amplitude  for  the  model  with  an  intermediate 
oceanic  path  of  31  km  length.  The  lower  plot  shows  the  corresponding  amplitude  for  the  model 
with  an  intermediate  oceanic  path  of  69  km  length.  The  horizontal  scale  between  the  two  vertical 
bars  is  different  for  each  plot.  This  region  represents  the  oceanic  path.  The  scale  between  the 
vertical  bars  is  chosen  so  that  the  shorter  oceanic  path  plots  with  the  same  length  as  the  longer 
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integral  will  be  discussed  in  the  next  section.  The  determination  of  displacement  and 
stress  Green’s  functions  is  necessary  in  all  applications  of  the  RT  integration  method 
regardless  of  whether  the  forcing  functions  are  FE  results,  analytic  stress  and  displace¬ 
ment  seismograms,  or  stresses  and  displacements  from  other  sources.  However,  the 

stress  seismograms  are  used  in  the  following  discussions  only  as  an  example  of  a  well 
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defined  form  of  forcing  functions.  Using  stress  and  displacement  seismograms  as  forc¬ 
ing  functions  produces  RT  integration  results  which  may  be  directly  compared  to 
purely  analytic  synthetics  allowing  one  to  verify  the  accuracy  of  the  RT  integration. 
For  the  SH  problem  in  Cartesian  coordinates  the  stresses  that  need  to  be  considered 
are  and  ertJ.  For  the  geometry  used  to  couple  surface  waves  from  a  FE  grid  into  a 
layered  medium  through  which  the  waves  will  be  transmitted  by  convolution  with 
propagator  matrix  generated  Green’s  functions,  only  the  stress  a is  used.  However, 
oXJ  will  also  be  derived  for  completeness.  Should  thfr  geometry  change  so  that  it 
would  be  necessary  to  integrate  over  a  horizontal  surface  such  as  the  bottom  of  the 
grid  then  aty  would  also  be  used. 

The  stresses,  <7^  and  oiy,  can  be  expressed  in  terms  of  spatial  derivatives  of  dis¬ 
placements.  One  method  of  calculating  values  for  these  stresses  is  to  express  the 
derivative  as  a  difference  and  evaluate  the  difference  equation  numerically.  The 
expressions  used  to  evaluate  the  stresses  in  the  FE  method  are  one  such  set  of 
difference  equations.  They  are 
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In  these  expressions  1,  2,  3,  4,  denote  the  positions  of  the  nodes  at  the  corners  of  the 
2-D  element  (Figure  2,  chapter  1)  for  which  the  element  center  stresses  are  <7xy  and 
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<r,r  The  displacement  perpendicular  to  the  FE  grid  is  uy,  and  the  horizontal  and  ver¬ 
tical  spacings  of  the  nodes  in  the  grid  are  Ax,  and  Az  respectively.  Evaluating  these 
relations  gives  a  reasonable  approximation  to  the  desired  stress  values.  A  second,  and 
more  direct,  method  for  determining  the  values  of  the  stresses  is  to  evaluate  the  ana¬ 
lytic  expression  for  each  stress  derived  from  the  analytic  displacement  expressions. 
The  derivations  of  analytic  stress  expressions  are  outlined  below.  Expressions  for  the 
stresses  are  determined  for  several  different  cases.  First,  dip  slip  and  strike  slip  double 
couple  sources  are  considered,  then,  a  line  source  in  a  horizontally  layered  medium. 
Finally,  the  expressions  for  the  line  source  Green’s  functions  are  derived. 

The  choice  of  evaluation  of  analytic  expressions  over  the  numerical  calculation  of 
derivatives  is  based  on  speed  and  accuracy.  Evaluation  of  the  analytic  expressions  for 
a  stress  component  requires  the  same  amount  of  calculation  as  the  evaluation  of  a  dis¬ 
placement.  For  a  given  depth  section  of  n  elements  n  element  center  stress  com¬ 
ponents  and  n  element  center  displacements  need  to  be  determined  to  evaluate  the  RT 
integral  which  propagates  the  wavefield  from  the  RT  surface  through  a  layered  struc¬ 
ture  to  the  receiver.  Thus,  if  the  calculation  of  each  stress  or  displacement  takes  time 
t,  the  computation  of  the  necessary  displacement  and  stress  components  by  evaluating 
the  analytic  expressions  would  take  time  2nt.  When  the  stress  component  is  deter¬ 
mined  using  numerical  differentiation,  time  (2n+2)t  is  required  to  evaluate  the  dis¬ 
placements  used  in  the  difference  equations,  time  nt  is  required  to  evaluate  the  element 
center  displacements,  and  additional  time  is  required  to  process  the  numerical 
difference  equations.  Clearly  direct  evaluation  is  faster.  The  calculation  of  numerical 
derivatives  is  known  to  be  a  potentially  unstable  numerical  procedure  since  subtrac¬ 
tion  of  almost  equal  numbers  is  possible.  Thus,  direct  evaluation  will  be  more  reliable 
in  cases  where  the  numerical  derivative  becomes  unstable.  The  analytic  evaluation  is 
also  more  accurate  since  it  is  equivalent  to  a  numerical  derivative  with  infinitesimal 


spacing  between  the  corners  of  the  element.  The  accuracy  of  the  numerical  derivative 
increases  as  that  spacing  decreases.  In  practice,  however,  for  the  types  of  disturbances 
considered  in  this  study,  the  increase  in  accuracy  is  small  and  both  estimates  are 
equally  acceptable. 

A  numerical  procedure  is  implemented  to  evaluate  the  analytic  expressions  for 
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the  stress  components.  The  validity  of  the  procedure  is  illustrated  by  comparing  ana¬ 
lytic  stress  time  histories,  determined  by  evaluation  of  equation  (4)  and  (6)  below,  with 
corresponding  results  generated  using  the  difference  «equations  (I).  The  results  of  these 
comparisons  are  discussed  after  the  analytic  stress  relations  are  developed. 

The  analytic  expressions  for  the  stress  components  for  SH  waves  from  a  point 
double  couple  source  are  directly  obtainable  from  the  expressions  used  to  determine 
the  corresponding  displacement  seismograms.  These  expressions, 

where  the  variables  are  as  defined  in  chapter  1,  lead  to  expressions  for  the  stress  <riy 
when  their  derivatives  with  respect  to  z  are  calculated.  Only  the  final  term  in  each 
equation  depends  directly  on  z.  From  Harkrider  (1964) 


d_  *r(z) 
dz  v0 


I  v0  /  cL  ]H 


Substituting  jv(z)J  from  (2)  for  u,  in  (1)  and  replacing  the  derivative  of 
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Expressions  (2)  also  yield  the  expressions  for  stress  a ^  when  their  derivatives 
with  respect  to  x  are  calculated.  All  terms  except  the  Hankel  function  are  constant 
with  respect  to  x.  By  expanding  the  Hankel  function  term  in  an  asymptotic  series  for 


large  r,  and  ignoring  terms  of  order  — ,  it  can  be  shown  that 
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Then  combining  (1),  (2),  and  (5)  gives 
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Equations  (2)  are  the  analytic  expressions  for  displacement  for  the  dip  slip  and  strike 
slip  faults,  the  corresponding  expressions  for  the  stresses  are  shown  in  (4)  and  (6). 


The  analytic  expressions  for  a  line  source  in  a  layered  medium  are  derived  by  a 
procedure  similar  to  that  used  in  the  previous  section  to  obtain  the  stress  expressions 
for  the  point  double  couple  source.  The  same  sequence  of  calculations  is  applied  to  the 
expression  for  the  displacement  produced  by  a  2-D  line  source  as  was  applied  to  the 
expressions  for  the  strike  slip  and  dip  slip  displacements.  These  calculations  yield  the 
analytic  2-D  line  source  stress  equations.  The  displacement  at  depth  i  due  to  a  line 
source  at  depth  h  is 
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Id  this  case  all  terms  in  (7)  except 


are  independent  of  2.  Thus,  the  stress 
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crtJ  has  the  same  form  as  the  displacement  expression  with 
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replaced  by  its 
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derivative  with  respect  to  z,  that  is  by  the  right  hand  expression  of  (3).  Similarly,  the 
stress  cTjy  has  the  same  form  as  the  displacement  with  the  exponential  in  x  replaced  by 
its  derivative.  Therefore,  the  stresses  for  the  2-D  line  source  are 
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Next  it  is  simple  to  extend  the  treatment  used  in  the  previous  paragraph  to  the 
expression  for  the  displacement  Green’s  function  for  a  line  source  in  a  layered  half 
space.  It  has  been  previously  shown  that  this  Green’s  function  is 
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In  this  case  a  stress  source  term  rather  than  a  stress  receiver  term  is  needed.  Thus, 
the  z  derivative  is  taken  with  respect  to  the  source  term.  All  terms  in  (9)  except 
vs0») 
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are  independent  of  z,  so  the  o„  term  has  the  same  form  as  the  displace- 
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ment  equation  with 
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replaced  by  its  derivative  with  respect  to  z.  The  form 
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of  this  derivative  is  identical  to  that  in  (3)  except  that  the  R  subscript  is  replaced  by 
an  S  subscript  denoting  the  properties  at  the  source  at  the  RT  surface  rather  than  at 
the  receiver.  Similarly  the  stress  <7X,  has  the  same  form  as  the  displacement  equation 
with  the  exponential  term  replaced  by  its  derivative  with  respect  to  x.  Therefore  the 
derivatives  of  the  Green’s  function  are 
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Now  it  is  necessary  to  illustrate  that  expressions  (4),  (6),  and  (8),  for  the  stress 
components  give  the  correct  results.  The  tests  illustrated  in  Figures  9,  10,  and  11  are 
analogous.  In  each  case  the  upper  trace  in  each  pair  is  the  displacement  or  stress  as 
calculated  using  one  of  the  sets  of  analytic  expressions  given  above.  The  lower  trace 
in  each  pair  is  the  same  displacement  component,  generated  by  averaging  the  sur¬ 
rounding  four  nodes,  or  the  same  stress  component  calculated  using  the  difference 
equations  (1)  with  Ax  *  Az  =■*  .5km.  The  node  spacings  are  chosen  to  correspond  to 
those  used  in  the  FE  calculations.  The  top  pairs  of  traces  are  the  V=uy  displace¬ 
ments  at  the  surface.  The  second  pairs  of  traces  are  the  stress  at  0.25  km  depth, 
and  the  third  pairs  of  traces  are  <riy  at  0.25  km  depth.  The  ratio  of  the  analytic  peak 
to  peak  amplitude  to  the  numerical  peak  to  peak  amplitude  is  given  for  each  pair  of 
traces  by  the  upper  number  beside  that  pair.  The  same  type  of  ratio,  analytic  to 
numerical,  is  calculated  using  the  RMS  amplitudes  of  the  seismograms  and  is  shown  as 
the  lower  number  beside  each  pair.  The  quality  of  the  waveform  correspondence  and 
the  amplitude  ratios  remain  essentially  constant  as  one  moves  down  the  depth  section. 
Figure  9  shows  the  results  for  the  case  of  the  2-D  line  source,  equation  (8).  The  pairs 
of  displacement  and  stress  time  histories  show  excellent  correspondence  when  their 
waveforms  are  compared.  Numerical  derivatives  for  agree  with  the  corresponding 
analytic  derivatives  to  within  %  4%  for  peak  to  peak  amplitude  and  to  within  «sl% 
for  RMS  amplitude.  The  peak  to  peak  differences  between  the  two  methods  are  about 
2%  for  the  average  displacement  V  and  the  stress  atJ.  As  expected,  for  all  three  cases, 
the  correspondence  improves  as  Ax  and  Ay  are  reduced.  Figure  10  shows  the  results 
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Figure  0:  Comparison  of  element  center  displacement  and  stress  time  histories 
calculated  for  a  2-D  line  source  using  two  methods.  Upper  traces  in  each  pair  are 
evaluations  of  the  analytic  expressions  at  the  element  centers.  Lower  traces  are  deter¬ 
mined  by  using  displacements  calculated  at  the  nodes  surrounding  the  element  center 
to  evaluate  the  difference  equations.  For  all  traces  A=sl500  km,  and  source  depth  is  8 
km.  The  fundamental  and  first  five  higher  modes  are  included  in  the  synthetics.  The 
upper  number  to  the  left  of  each  trace  is  the  peak  to  peak  amplitude  ratio  of  the 
numerical  to  analytic  calculations,  the  lower  number  is  the  corresponding  RMS  ampli¬ 
tude  ratio. 
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for  the  case  of  a  strike  slip  source,  equations  (4)  and  (6),  using  the  fundamental  mode 
only.  Again,  the  waveform  correspondence  within  each  pair  is  excellent.  The  peak  to 
peak  amplitude  ratios  are  close  to  one,  showing  less  than  4%  difference.  Since  the 
major  variation  between  the  two  methods  is  the  resolution  of  the  magnitude  of  the 
high  frequency  peak  near  the  onset  of  the  trace  the  RMS  amplitudes  show  a  smaller' 
variation  of  2.5%  or  less.  Figure  11  shows  the  results  for  the  strike  slip  source,  for  a 
mode  sum  including  the  fundmental  mode  and  the  first  five  higher  modes.  The 
waveforms  do  not  match  as  precisely  as  in  the  case  of  the  fundamental  mode  alone, 
but  the  correspondence  is  still  excellent.  The  amplitude  ratios  do  not  significantly 
change  when  the  higher  mode  energy  is  added.  Again,  both  for  the  fundamental  mode 
alone  and  for  the  mode  sum,  correspondence  between  the  analytic  and  numerical 
methods  of  evaluation  improves  rapidly  as  Ax  and  Az  are  reduced.  These  tests  not 
only  indicate  the  validity  of  the  analytic  expressions  above,  they  also  give  an  estimate 
of  the  uncertainties  present  in  the  numerical  differentiation  used  in  the  FE  calcula¬ 
tions. 

RT  Coupling  of  Analytic  2-D  Seismograms  and  Green’s  Functions 
In  this  section  the  validity  and  accuracy  of  the  numerical  implementation  of  the 
Representation  Theorem  coupling  technique,  used  in  the  following  section  to  pass  FE 
results  into  a  layered  media,  will  be  discussed  in  detail.  The  tests  discussed  below  are 
designed  to  give  RT  integration  results  directly  comparable  with  purely  analytic 
results.  First  mode  by  mode  results  are  presented  to  illustrated  where  the  discrepan¬ 
cies  between  the  RT  results  and  the  analytic  results  originate.  Mode  sum  results  are 
then  presented.  The  form  of  the  RT  integral  and  its  relation  to  the  propagator  matrix 
formulation  is  discussed  and  a  method  for  mode  by  mode  filtering  of  mode  sum  forcing 
function  input  is  explained.  The  type  of  formulation  used  to  explain  the  filtering 
method  is  also  applied  to  derive  more  quantitative  mode  by  mode  analysis  of  the 
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origins  of  discrepancies  between  RT  integration  results  and  purely  analytic  solutions 
These  quantitative  estimates  of  the  sources  of  discrepancies  explain  the  differences 
seen  between  the  RT  integration  results  and  the  purely  analytic  results  quite  well. 

All  the  results  discussed  in  this  section  are  derived  using  a  simple  geometry.  In 

% 

all  cases  the  model  is  a  layer  over  a  half-space.  The  layer  has  a  thickness  of  thirty 
two  kilometers,  a  SH  wave  velocity  3.5  km/s,  and  a  density  2.7  g/cm3  The  half-space 
has  SH  wave  velocity  4.5  km/s  and  density  3.4  g/cm3.  The  same  layer  over  a  half 
space  model  is  used  for  the  entire  path,  making  calculation  of  purely  analytic  synthetic 
for  the  entire  path  length  of  A  =  1600  km,  or  A  =  1750  km  simple.  Purely  analytic 
synthetic  seismograms  were  calculated  at  these  distances  for  the  mode  sum  and 
separately  for  each  of  the  fundamental  mode  and  the  first  five  higher  modes.  The 
forcing  functions  used  are  the  displacement  and  stress  seismograms  for  a  line  source  at 
at  depth  of  ten  kilometers  and  a  distance  of  A—  1500  km  from  the  source.  The  forc¬ 
ing  functions  are  evaluated  at  positions  corresponding  to  the  element  centers  of  the 
rightmost  column  of  elements  in  a  FE  grid  with  horizontal  and  vertical  spacing  of  5 
km,  whose  right  hand  edge  lies  A  —  1500.25  km  from  the  source.  Thus,  the  seismo¬ 
grams  are  evaluated  at  points  along  a  vertical  surface,  at  depth  intervals  of  0.5  km, 
beginning  at  a  depth  of  0.25  km  below  the  surface  All  forcing  functions  are  evaluated 
at  a  distance  of  A  =  1500  km.  Separate  sets  of  forcing  functions  were  generated  for 
each  mode.  A  set  of  forcing  functions  which  is  a  sum  over  the  fundamental  and  the 
first  five  higher  modes  was  also  calculated.  The  displacement  seismograms  calculated 
using  the  RT  integral  are  for  distances  of  A  =  1600  km,  orA  =  1750  km.  Performing 
RT  integrations  to  determine  hybrid  solutions  at  these  distances  requires  Green’s  func¬ 
tions  for  propagation  distances  A  —  100,  or  A  =  250  km.  At  each  of  these  distances 
a  set  of  stress  and  displacement  Green’s  functions  was  determined  for  each  mode  and 
an  additional  set  was  determined  for  a  mode  sum  including  the  fundamental  and  the 
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first  five  higher  modes.  The  Green’s  functions  in  each  set  were  evaluated  for  a  source 
at  each  of  the  locations  where  displacement  and  stress  forcing  functions  were  deter* 
mined  and  a  receiver  at  the  surface.  The  RT  integration  surface  for  these  examples 
extended  to  a  depth  of  37.5  km  and  included  seventy  five  integration  points.  In  the 
following  discussions  the  seismogram  resulting  from  a  RT  integration  may  be  referred 
to  as  a  hybrid  seismogram,  even  though,  for  these  testa  of  accuracy,  the  same  method 
is  being  used  to  generate  Green’s  functions  and  forcing  functions. 

The  first  group  of  tests  using  the  sets  of  forcing  functions  and  Green’s  functions 
discussed  above  produced  mode  by  mode  RT  integration  results  to  compare  to  the 
purely  analytic  synthetic  single  mode  seismograms.  The  set  of  single  mode  forcing 
functions  for  each  of  the  individual  modes,  was  convolved  with  the  single  mode  set  of 
Green’s  functions  for  the  same  mode  according  to  the  RT  integration  relation.  This 
produced  a  hybrid  seismogram  for  that  mode  to  compare  to  the  corresponding  purely 
analytic  synthetic.  Comparisons  of  the  RT  integration  sums  and  the  purely  analytic 
synthetics  for  each  individual  mode  are  shown  in  Figures  11  and  12.  In  Figure  11 
comparisons  for  the  fundamental  mode,  the  first  higher  mode,  and  the  second  higher 
mode  are  presented.  Each  of  these  pairs  is  illustrated  at  a  distance  of  A  =  1600  km 
from  the  source.  Figure  12  is  a  continuation  of  Figure  11  showing  the  results  for  the 
third  through  fifth  higher  modes.  The  upper  two  pairs  of  traces  are  also  illustrated  for 
a  distance  A  =  1600  km  .  The  lowermost  pair  of  traces,  those  for  the  fifth  higher 
mode,  are  illustrated  at  a  distance  of  A  =»  1750  km  from  the  source.  The  fifth  higher 
mode  results  at  1500  km  are  equally  well  fit,  but  presenting  an  example  for  another 
distance  helps  verify  the  observation  that  the  goodness  of  fit  between  the  hybrid 
seismograms  and  the  analytic  synthetics  does  not  depend  significantly  on  distance  pro¬ 
pagated  using  the  RT  integration  over  stress  and  displacement  seismograms  and 
Green’s  functions.  In  each  of  these  figures  there  are  three  pairs  of  seismograms.  The 
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Figure  12:  Comparison  of  analytic  and  hybrid  synthetics  at  A=*1600km.  The 
hybrid  synthetics  are  transmitted  the  first  1500  km  using  the  propagator  method,  and 
the  remaining  100  km  using  RT  integration  with  propagator  generated  Green’s  func¬ 
tions.  The  left  column  of  seismograms  shows  the  analytic  synthetics  and  the  right 
column  the  hybrid  synthetics.  The  results  are  presented  mode  by  mode,  and  the  mode 
is  identified  below  each  pair  of  traces.  The  number  between  each  pair  of  traces  is  the 
ratio  of  the  RMS  amplitude  of  the  hybrid  trace  to  the  RMS  amplitude  of  the  analytic 
trace. 
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Figure  13:  Comparison  of  analytic  and  hybrid  synthetics  at  £=1600  km.  This 
figure  is  a  continuation  of  Figure  12  showing  the  remaining  modes.  Other  details  are 
the  same  as  for  Figure  12. 
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leftmost  seismogram  in  each  pair  shows  the  purely  analytic  synthetic,  the  rightmost 
seismogram  shows  the  hybrid  synthetic  produced  by  RT  integration.  Between  each 
pair  of  seismograms  is  a  number  indicating  the  ratio  of  the  RMS  amplitude  of  the 
hybrid  synthetic  to  the  RMS  amplitude  of  the  purely  analytic  synthetic.  Below  each 
pair  of  traces  is  a  label  indicating  which  mode  is  being  illustrated.  All  the  seismo¬ 
grams  are  bandpass  filtered  for  periods  between  one  and  twenty-five  seconds.  The 
long  period  limit  on  the  band  pass  filter  was  chosen  to  improve  the  correspondence 
between  the  waveform  of  the  purely  analytic  synthetic  and  the  RT  integration  result 
for  the  same  mode.  Longer  periods  seem  to  be  poorly  reconstructed  by  the  RT 
integration,  and  are  consequently  filtered  out  of  the  displayed  results.  The  filtering 
has  the  largest  effect  for  the  first  two  higher  modes,  and  has  a  progressively  smaller 
effect  for  each  successive  higher  mode.  Higher  mode  results  depend  less  on  the  long 
periods  and  produce  an  increasingly  good  fit  even  before  filtering.  The  increasing 
discrepancies  between  the  purely  analytic  and  the  hybrid  synthetics  for  successively 
lower  higher  modes  is  due  to  the  larger  proportion  of  long  periods  in  those  modes,  and 
is  a  likely  source  of  discrepancies  in  the  mode  sum  calculations  presented  in  the  next 
paragraph  and  will  be  explained  later  in  terms  of  the  quantitative  estimates  of  error 
yet  to  be  derived. 

The  next  test  conducted  was  the  RT  integration  using  the  mode  sum  forcing 
functions  and  Greens  functions.  The  hybrid  synthetic  resulting  from  this  integration 
is  compared  to  the  purely  analytic  mode  sum  synthetic  in  Figure  14.  In  this  figure  the 
hybrid  synthetic  is  labeled  mode  sum  and  the  purely  analytic  synthetic  is  labeled  ana¬ 
lytic.  The  seismogram  labeled  sum  of  single  modes  is  calculated  by  summing  the 
hybrid  solutions  for  each  individual  mode  to  produce  another  estimate  of  the  mode 
sum  hybrid  synthetic.  Clearly  the  waveforms  of  all  three  seismograms  are  extremely 
similar.  The  numbers  to  the  left  of  each  of  the  lower  two  seismograms  indicates  the 
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ratio  of  the  RMS  amplitude  of  that  seismogram  to  the  RMS  amplitude  of  the  purely 
analytic  mode  sum  seismogram.  The  agreement  in  amplitude  between  the  purely  ana¬ 
lytic  and  hybrid  mode  sum  results  is  better  than  the  agreement  seen  for  any  single 
mode.  This  improvement  in  agreement  is  probably  fortuitous.  The  mode  sum  syn¬ 
thetic  calculated  as  a  sum  over  the  single  mode  hybrid  results  shows  a  more  realistic 
amplitude  measure  consistent  with  the  single  mode  results  previously  presented. 

The  RT  integration  is  accomplished  by  the  evaluation  of  the  following  expression 
based  on  integrating  equation  32  of  chapter  1  along  a  vertical  surface. 
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In  this  expression  the  displacement  at  time  t  at  location  7  is  determined  as  a  RT 
integral.  The  forcing  functions,  and  Psu2,i(fi>^3)  &re  evaluated  on  the  verti¬ 

cal  surface  perpendicular  to  the  propagation  direction  at  a  distance  from  the  origin 
for  points  with  a  range  of  values  of  £3.  Thus,  the  propagation  distance  for  the  forcing 
functions  is  ft.  The  Green’s  functions,  r»  and  rM,i  are  evaluated  by  placing  a  source 
point  in  each  position  where  a  forcing  function  is  evaluated.  By  substituting  alternate 
expressions,  in  terms  of  the  variables  used  in  the  Propagator  matrix  method,  for  the 
Green’s  functions  and  forcing  functions  in  equation  (11)  the  Representation  Theorem 
can  be  expressed  in  terms  of  terms  constant  with  respect  to  (3  times  a  simple  integral 
with  £3  as  an  integration  variable.  To  derive  this  form  of  the  RT  equation  (9)  is  sub¬ 
stituted  for  r^x.z;^,^),  and  equation  (8b)  is  substituted  for  r22,i(x,z;ft,ft).  Also, 
equation  (7)  is  substituted  for  Uj(ft,ft),  and  equation  (10b)  for  u2>i(ft,ft).  Performing 
these  substitutions  gives 
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When  the  above  mentioned  substitutions  into  equation  (llj  ar  performed  it  becomes 
evident  that  both  convolutions  in  the  integrand  of  equation  (11)  produce  identical 
expressions.  Therefore,  it  should  be  possible  to  accelerate  the  numerical  evaluation  of 
*11)  by  evaluation  only  one  of  the  convolutions  and  then  doubling  the  resulting  solu- 
tion.  This  approach  had  the  additional  advantage  that  it  makes  it  unecessary  to 
record  both  displacement  and  stress  time  histories  in  the  FE  calculations.  Either  one 
of  these  should  be  sufficient  to  calculate  a  2-D  SH  RT  integral.  Even  more  useful  is 
the  fact  that  the  equality  of  the  two  convolution  terms  allows  determination  of  the 
value  of  the  RT  integral  with  the  evaluation  of  half  the  number  of  Greens  function. 
Tests  have  been  conducted  to  investigate  the  validity  of  this  approach.  The 
waveforms  of  the  solutions,  for  complete  evaluation  of  the  RT  integral,  and  for  evalua¬ 
tion  of  either  term  in  the  integral  are  indistinguishable  when  the  resulting  seismograms 
are  examined.  Thus,  it  was  considered  unnecessary  to  illustrate  the  results  from  these 
tests.  The  amplitudes  of  the  sums  of  each  term  in  the  integral  of  (11),  however,  were 
ot  necessarily  identical.  Using  the  term  containing  the  displacement  and  the  Green’s 
function  stress  general  yielded  amplitudes  a  percent  or  two  higher  than  using  the  other 
term.  These  amplitudes  were  usually  in  better  agreement  with  the  synthetic  than  the 
amplitudes  determined  by  summing  the  two  terms.  The  differences  in  amplitude 
agreement  are  small  enough  to  be  ignored.  Further  calculations  of  the  RT  integral 
may  be  done  by  doubling  the  value  of  the  first  term  in  (11)  without  significantly 
effecting  the  solution. 

Returning  to  equation  (12)  it  is  clear  that  many  of  the  terms  do  not  depend  upon 
the  integration  variable  f3.  Taking  these  terms  outside  the  integral  and  then  compar¬ 
ing  them  to  equation  (7)  allows  equation  (12)  to  be  simplified  to 
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immediately  leads  to  the  relation 
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Later  the  use  of  this  last  relation  in  the  estimation  of  the  minimum  error  for  each 
mode  will  be  discussed.  First  some  consequences  of  equation  (14)  will  be  considered. 
It  is  well  known  that  if  i  and  j  represent  two  different  modes  for  a  given  period,  that  is 
k^kj  for  ,  or  for  a  given  wave  number  kj*kj  ,  then  the  orthogonality 

relation  for  Love  waves  states  that 


f  p(z)vi(z)Yj(z)dz  ^  0  iy*j  (16) 

o 

In  this  equation  Vj  is  the  equivalent  of  the  component  of  u3  in  equation  (14)  due  to  the 
single  mode  i.  Comparing  (14)  and  (16)  shows  that  equation  (14)  is  a  form  of  the 


orthogonality  relation.  At  this  point  it  is  useful  to  notice  that  the  two 
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terms  in  the  equation  (14)  each  have  separate  origins.  One  originates  with  the  forcing 
functions  and  the  other  with  the  Green's  functions.  Thus,  any  single  modes  not  com¬ 
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from  the  Green’s  function  will  produce  zero  contribution  to  the  resulting  hybrid  result. 
This  implies  that  the  only  modes  present  in  both  the  Green’s  functions  and  the  forcing 
functions  will  be  present  in  the  RT  integration  results.  Thus,  choosing  Green’s  func¬ 
tions  with  a  subset  of  the  modes  present  in  the  forcing  function  will  produce  a  filter 
that  gives  RT  integration  results  that  contain  only  that  subset  of  modes.  This  is 


demonstrated  in  the  next  set  of  numerical  experiments. 

A  series  of  calculations  investigating  the  accuracy  and  efficiency  of  the  use  of  a 
Green’s  functions,  containing  only  a  subset  of  the  modes  present  in  the  forcing  func¬ 
tions,  as  a  filter  to  extract  only  those  modes  from  the  forcing  functions  has  been  com¬ 
pleted.  In  particular,  the  single  mode  sets  of  Green's  functions  were  used  in  the  RT 
integral  along  with  the  mode  sum  set  of  forcing  functions.  A  RT  integration  was  com¬ 
pleted  for  each  single  mode  set  of  Green’s  functions.  Figure  15  presents  results  of  the 
RT  integration  of  the  mode  sum  forcing  functions  and  the  fundamental  mode  Green’s 
functions.  Figure  16  is  analogous  to  figure  15  but  uses  the  third  higher  mode  Green’s 
function  set.  In  each  of  these  figures  two  columns  of  five  seismograms  are  illustrated. 
The  left  column  shows  results  of  the  single  mode  Green’s  functions  integrated  with  the 
mode  sum  forcing  functions.  The  right  column  shows  the  results  of  integrating  single 
mode  Green’s  functions  with  single  mode  forcing  functions.  The  RT  integration  is 
numerically  evaluated  by  summing  the  convolutions  at  discrete  points  along  the  RT 
surface.  The  sum  begins  with  the  point  nearest  the  surface  and  then  adds  points  at 
steadily  increasing  depth.  In  Figures  15  and  16  the  first  row  of  seismograms  is  a  single 
convolution,  the  sum  down  to  a  depth  of  0.25  km.  The  second  row  is  the  sum  of  15 
convolutions,  and  includes  all  integration  points  to  a  depth  of  7.25  km.  This  pattern 
continues  with  the  depth  of  the  deepest  point  included  in  the  integration  indicated  to 
the  right  of  each  pair  of  seismograms.  The  numbers  between  each  pair  of  seismograms 
indicate  the  ratio  of  the  RMS  amplitude  of  the  left  trace  to  the  RMS  amplitude  of  the 
right  trace.  For  the  fundamental  mode  case  this  amplitude  ratio  varies  only  slightly 
with  depth  and  shows  a  trend  of  increasing  with  depth  only  in  the  upper  eight  kilome¬ 
ters.  However,  as  illustrated  by  the  third  higher  mode,  the  higher  mode  ratios  con¬ 
tinue  increasing  with  depth  to  much  larger  depths.  This  implies  that  the  integration 
must  proceed  to  a  reasonable  depth  for  the  amplitude  of  the  filtered  mode  sum 


Figure  15:  Comparison  of  hybrid  synthetics  at  A— 1600  km.  The  traces  in  the 
left  column  are  calculated  by  ET  convolution  of  the  mode  sum  forcing  functions  calcu¬ 
lated  at  A— 1500  km  with  the  fundamental  mode  Green’s  functions  for  100  km 
further  propagation.  The  traces  in  the  right  column  are  determined  by  RT  convolu¬ 
tion  of  fundamental  mode  forcing  functions  for  DELTA-1500  km  and  the  same 
Green’s  functions  as  the  left  column.  The  numbers  to  the  right  of  each  pair  of  seismo¬ 
grams  show  the  depth  to  which  integration  along  the  RT  integration  surface  has  been 
completed.  The  central  numbers  show  the  RMS  amplitude  ratio  of  the  the  right  trace 
to  the  left  trace. 
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Figure  18:  Comparison  of  hybrid  synthetics  at  A >—1600  km.  The  traces  in  the 
left  column  are  calculated  by  RT  convolution  of  the  mode  sum  forcing  functions  calcu- 
lated  at  A=— 1500  km  with  the  third  higher  mode  Green’s  functions  for  100  km  further 
propagation.  The  traces  in  the  right  column  are  determined  by  RT  convolution  of 
third  higher  mode  forcing  functions  with  the  same  Green’s  functions  as  the  left 
column.  The  numbers  to  the  right  of  each  pair  of  seismograms  and  between  each  pair 
of  seismograms  have  the  same  meanings  as  for  Figure  14. 
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seismogram  to  reach  that  of  the  single  mode  hybrid  results.  For  the  modes  investi¬ 
gated  it  was  determined  that  the  amplitude  ratios  stabilized  near  or  above  the  base  of 
the  crustal  layer,  and  that  when  the  amplitudes  had  stabilized  the  waveforms  of  the 
two  calculations  showed  little  difference  and  remained  relatively  stable.  Figures  17 
and  18  illustrate  the  mode  by  mode  results  of  the  test  discussed  above.  Rather  than 
show  the  detailed  progression  to  the  stable  result  as  in  the  previous  two  figures,  a 
depth  of  thirty  kilometers  was  chosen  as  the  vertical  extent  of  integration.  This 
allows  the  illustration  of  the  results  of  all  six  modes  in  two  figures.  In  these  figures 
two  columns  of  seismograms  are  illustrated.  The  left  column  shows  the  results  using 
the  mode  sum  forcing  functions  and  the  single  mode  Green’s  functions.  The  right 
column  shows  the  results  using  the  single  mode  forcing  functions  and  the  single  mode 
Green’s  functions.  The  numbers  between  each  pair  of  seismograms  indicate  the  ratio 
of  the  RMS  amplitude  of  the  left  trace  to  the  RMS  amplitude  of  the  right  trace 
Below  each  pair  of  seismograms  the  mode  of  the  Green’s  functions  is  indicated.  These 
seismograms  appear  different  from  those  shown  in  Figures  12  and  13  since  they  have 
not  been  filtered  before  plotting.  In  summary,  the  effect  of  using  a  set  of  Green's 
functions  containing  a  subset  of  the  modes  present  in  the  set  of  forcing  functions  is  to 
produce  an  efficient  filter  that  allows  only  the  modes  common  to  both  sets  to  appear 
in  the  hybrid  result. 

Now  it  is  time  to  return  to  the  analysis  of  the  accuracy  of  the  RT  integration  for 
a  given  single  frequency  mode.  Rearranging  equation  (15)  gives 

2AtIt-l  (I7) 

Evaluation  of  this  simple  equation  provides  a  direct  estimate  of  the  accuracy  of  the 
integration  on  a  mode  by  mode  basis.  The  estimate  of  the  accuracy  is  obtained  by 
evaluating  the  I,  integrand  at  each  integration  point  used  in  the  RT  integration,  for 
each  frequency  on  each  branch  of  the  dispersion  curve  used.  For  each  single  frequency 
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Figure  17;  Comparison  of  synthetics  at  A»»1600km.  The  synthetics  are 
transmitted  the  first  1500  km  using  the  propagator  method,  and  the  remaining  100  km 
using  RT  integration  with  propagator  generated  Green’s  functions.  The  left  column  of 
seismograms  shows  the  RT  convolution  of  the  mode  sum  forcing  functions  and  the 
Green’s  functions  of  the  indicated  mode,  the  right  column  shows  synthetics  resulting 
from  RT  integration  of  forcing  functions  and  Green’s  functions  of  the  indicated  mode 
only.  The  results  are  presented  mode  by  mode,  and  the  mode  is  identified  below  each 
pair  of  traces.  The  number  between  each  pair  of  traces  is  the  ratio  of  the  RMS  ampli- 
tude  of  the  left  trace  to  right  trace. 
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Figure  18:  Comparison  of  synthetics  at  A»1600  km.  This  figure  is  a  continua¬ 
tion  of  Figure  12  showing  the  remaining  modes.  Other  details  are  the  same  as  for  Fig¬ 
ure  12. 


mode  the  quantity  on  the  right  hand  side  of  equation  17  is  determined  at  each  integra¬ 
tion  point  and  summed  over  the  integration  surface.  If  the  solution  were  perfect  with 
no  error  present,  then  the  sum  would  be  exactly  one.  In  practice  the  sum  departs  from 
one  by  some  amount  which  gives  and  estimate  of  the  size  of  the  minimum  error  that 
could  be  expected  in  that  mode  in  the  RT  integration  results.  The  estimate  is  a,, 
minimum  since  it  does  not  account  for  the  phase  of  the  arrivals  nor  for  possible  errors 
in  that  phase.  The  evaluation  of  the  error  using  this  relation  is  much  faster  than  com¬ 
paring  results  from  multiple  applications  of  RT  theorem  coupling.  For  example  gen¬ 
erating  the  tables  on  pages  219  and  220  which  show  the  effects  of  varying  grid  spacing 
or  of  varying  the  vertical  extent  of  tbe  RT  integration  surface  on  the  errors  took  less 
than  25%  the  time  needed  for  a  single  RT  integration  example. 

The  tables  on  pages  219  and  220  show  the  effects  of  grid  spacing  and  vertical 
extent  of  the  integration  surface.  In  the  tables  z  indicates  the  depth  in  kilometers  to 
which  the  integration  surface  extends,  Ax  indicated  the  spacing  between  successive 
integration  points  on  that  surface,  and  T  indicates  the  period  of  the  given  mode.  For 
each  combination  of  these  parameters  illustrated  the  numbers  in  the  table  are  the 
values  of  2Ai.Ii-  It  is  clear  from  examining  these  tables  that  the  grid  spacing  is  fine 
enough  to  yield  good  results.  Halving  or  doubling  the  grid  spacing  makes  less  than 
one  percent  difference  in  the  error  estimates  which  are  at  a  level  of  about  98.5%. 
Reducing  the  time  spacing  by  as  much  as  a  factor  of  four  made  no  significant 
differences  in  the  estimates  of  accuracy.  So  no  table  was  included  to  illustrate  how 
accuracy  changed  with  time  spacing.  Seeing  the  errors  on  a  frequency  by  frequency 
basis  also  helps  to  understand  the  causes  of  the  inaccuracies  noted  in  the  actual  RT 
integration  results.  The  only  variation  that  made  noticeable  changes  in  tbe  accuracy 
of  individual  modes  was  variation  in  the  depth  to  which  tbe  integration  surface  was 
extended.  This  change  was,  as  expected,  largest  in  the  longest  period  modes.  At 
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0.9851 

0.9851 
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0.9851 
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0.9850 
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longer  periods  the  wavelength  increases  and  the  integration  surface  needs  to  be 
extended  downward  to  maintain  accurate  results.  When  the  depth  to  which  the 
integration  surface  extends  is  less  than  n  wavelengths  then  significant  error  occurs. 
This  relation  is  clearly  visible  in  the  tables  on  pages  222  and  223.  These  tables  show 
the  mode  by  mode  values  of  2AiJi  for  each  separate  overtone,  and  for  the  fundament 
tal.  It  is  easily  seen  that  the  error  is  very  small  until  n+1  wavelengths,  for  the  nth 
higher  mode,  become  longer  than  the  depth  extent  of  the  RT  integration  surface. 
This  explains  why  the  RT  integration  results  and  the  analytic  synthetics  needed  to 
have  the  long  periods  removed  to  yield  a  good  fit.  Those  long  period  components  had 
large  errors  due  to  incomplete  representation  caused  by  the  truncation  of  the  integra¬ 
tion  surface  at  depth  and  thus  contaminated  otherwise  accurate  results. 

Using  the  RT  Integration  Method  for  Long  Oceanic  Paths 

Now  that  the  reliability  of  the  RT  coupling  of  FE  results  to  a  receiver  distant 
from  the  FE  grid  has  been  demonstrated  it  can  be  applied  to  the  problem  of  investi¬ 
gating  the  propagation  of  Lg  when  the  intermediate  oceanic  path  is  too  long  to  be 
handled  using  FE  alone.  The  amount  of  calculation  needed  to  determine  each  seismo¬ 
gram  using  the  coupling  technique  is  large  enough  that  the  calculation  of  an  entire 
depth  section  for  reentry  into  another  FE  calculation  would  be  worthwhile  only  if  a 
specific  example  were  to  be  considered,  and  detailed  structure  were  known.  Thus,  the 
approach  taken  here  is  to  calculate  seismograms  at  the  surface  at  receivers  distant 
from  the  end  of  the  FE  grid,  and  to  examine  what  the  remaining  modal  energy  looks 
like  after  propagation  for  a  large  distance  through  the  oceanic  crust.  Then  by  know¬ 
ing  the  proportion  of  surface  amplitude  that  is  removed  by  transmission  through  a 
reverse  transition  an  estimate  of  the  RMS  amplitude  of  the  Lg  energy  that  travels 
through  the  transition  can  be  made. 
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14.00 

0.9701 

8.400 

0.8488 

4.200 

0.9748 

12.00 

0.9782 

8.200 

0.8729 

4.100 

0.9778 

10.00 

0.9827 

8.000 

0.8935 

4.000 

0.9800 

8.000 

0.9848 

7.800 

0.9108 

3.800 

0.9827 

6.000 

0.9854 

7.600 

0.9254 

3.600 

0.9842 

4.000 

0.9852 

7.400 

0.9374 

3.400 

0.9849 

2.000 

0.9848 

7.200 

0.9474 

3.200 

0.9852 

1.800 

0.9848 

7.000 

0.9555 

3.000 

0.9852 

1.600 

0.9847 

6.500 

0.9695 

2.500 

0.9851 

1.400 

0.9847 

6.000 

0.9775 

2.000 

0.9849 

1.200 

0.9846 

5.500 

09818 

1.900 

0.9849 

1.000 

0.9847 

5.000 

0.9839 

1.800 

0.9849 

0.8000 

0.9847 

4.500 

0.9848 

1.700 

0.9848 

0.7000 

0.9847 

4000 

0.9851 

1.600 

0.9848 

0.6000 

0.9839 

3.500 

0.9852 

1.500 

0.9848 

0.5000 

0.9847 

3.000 

0.9851 

1.400 

0.9848 

0.4000 

0.9847 

2.500 

0.9850 

1.300 

0.9847 

0.3000 

0.9848 

2.000 

0.9848 

1.200 

0.9847 

0.2500 

0.9889 

1.800 

0.9848 

1.100 

0.9847 

0.2000 

0.9849 

1.600 

0.9847 

1.000 

0.9847 

0.1500 

0.9849 

1.400 

0.9847 

0.9001 

0.9846 

0.1000 

0.9851 

1.200 

0.9846 

0.8001 

0.9846 

0.0000 

0.0000 

1.000 

0.9846 

0.7001 

0.9846 

0.0000 

0.0000 

0.9001 

0.9846 

0.6001 

0  9844 

0.0000 

0.0000 

0.8001 

0.9845 

0.5001 

0.9843 

0.0000 

0.0000 

0.7001 

0.9849 

0.4501 

0.9844 

0.0000 

0.0000 

0.6001 

0.9843 

0.4000 

0.9845 

0.0000 

0.0000 

0.5001 

0.9844 

0.3500 

0.9845 

0.0000 

0.0000 

0.4001 

0.9841 

0.3000 

0.9845 

0.0000 

0.0000 

0.3001 

0.9841 

0.2500 

0.9849 

0.0000 

0.0000 

0.2501 

0.9844 

0.2000 

0.9850 

0.0000 

0.0000 

0.2000 

0.9849 

0.1500 

0.9851 

0.0000 

0.0000 

0.1500 

0.9851 

0.100E 

0.9865 

0  0000 

0.0000 

0.1001 

4.234 

0.0000 

0  0000 

3.831 

0.0000 

2.873 

0.0000 

2.299 

0.0000 

3.824 

.07248 

2.869 

.00512 

2.296 

0.1171 

3.808 

0.2003 

2.860 

0.2664 

2.290 

0.3185 

3.800 

0.2651 

2.854 

0.3519 

2.286 

0.4144 

3.700 

0.6640 

2.800 

0.7295 

2.200 

0.9187 

3.600 

0.8238 

2.700 

0.9145 

2.100 

0.9748 

3.500 

0.0005 

2.600 

0.9624 

2.000 

0.9839 

3.400 

0.0400 

2.500 

0.9777 

1.000 

0.9855 

3.300 

0.0600 

2.400 

0.9830 

1.800 

0.0857 

3.200 

0.0722 

2.300 

0.9848 

1.700 

0.0856 

3.100 

0.0784 

2.200 

0.9854 

1.600 

0.0854 

3.000 

0.0818 

2.100 

0.9855 

1.500 

0.9853 

2.000 

0.0836 

2.000 

0.9855 

1.400 

0.0852 

2.800 

0.0846 

1.000 

0.9854 

1.300 

0.0851 

2.700 

0.0851 

1.800 

0.9853 

1.200 

0.9850 

2.000 

0.8853 

1.700 

0.9852 

1.100 

0.9850 

2.500 

0.0854 

1.600 

0.9851 

1.000 

0.9849 

2.400 

0.08S4 

1.500 

0.9850 

0.9501 

0.9849 

2.300 

0.0853 

1.400 

0.0850 

0.9001 

0.9849 

2.200 

0.0853 

1.300 

0.9849 

0.8501 

0.9849 

2.000 

0.0851 

1.200 

0.9849 

0.8001 

0.9849 

1.800 

0.0850 

1.100 

0.9848 

0.7501 

0.9848 

1.600 

0.0840 

1.000 

0.9848 

0.7001 

0.9848 

1.500 

0.0840 

0.9501 

0.0848 

0.6501 

0.9848 

1.400 

0.0848 

0.9001 

0.9848 

0.6000 

0.9848 

1.300 

0.0848 

0.8501 

0.9848 

0.5500 

0.9848 

1.200 

0.0848 

0.8001 

0.0848 

0.5000 

0.9848 

1.100 

0.0848 

0.7501 

0.9847 

0.4500 

0.9848 

1.000 

0.0847 

0.7001 

0.0847 

0.4000 

0.9847 

0.0001 

0.0847 

0.6501 

0.9847 

0.3501 

0.9847 

0.8001 

0.0847 

0.6000 

0.9847 

0.3000 

0.9847 

0.7501 

0.0846 

0.5500 

0.9847 

0.2500 

0.9851 

0.7001 

0.0847 

0.5000 

0.9847 

0.2000 

0.9852 

0.6501 

0.0846 

0.4500 

0.9847 

0.1500 

0.9846 

0.6001 

0.0846 

0.4000 

0.9847 

0.1000 

0.9849 

0.5500 

0.9846 

0.3500 

0.9846 

0.0000 

0.0000 

0.5000 

0.9846 

0.3000 

0.9846 

0.0000 

0.0000 

0.4500 

0.9846 

0.2500 

0.9851 

0.0000 

0.0000 

0.4000 

0.9840 

0.2000 

0.9846 

0.0000 

0.0000 

0.3500 

0.9846 

0.1500 

0.9847 

0.0000 

0.0000 

0.3000 

0.9845 

0.1000 

0.9849 

0.0000 

0.0000 

0.2500 

0.9844 

0.0000 

0.0000 

0.0000 

0.0000 

0.2000 

0.9837 

0.0000 

0.0000 

0.0000 

0.0000 

0.1500 

0.9853 

0.0000 

0.0000 

0.0000 

0.0000 

0.1000 

0.9856 

0.0000 

0.0000 

0.0000 

0.0000 
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In  order  to  perform  the  RT  integration  that  propagates  the  FE  results  across  the 
oceanic  portion  of  the  path  sets  of  Green’s  functions  must  first  be  calculated  for  each 
distance  at  which  the  seismogram  is  to  be  examined.  It  should  be  noted  that  these 
Green’s  Functions  are  not  identical  to  those  used  in  the  tests  above.  Since  the  propa¬ 
gation  path  has  an  oceanic  crustal  structure  a  different  set  of  modes  appropriate  for  an* 
oceanic  structure  rather  than  a  continental  one  must  be  used  to  determine  the  Green’s 
functions.  Since  the  oceanic  crust  is  thinner  than  the  continental  crust  the  modes 
tend  to  concentrate  at  higher  frequencies.  Thus,  more  (eight)  higher  modes  were  used 
in  the  oceanic  mode  set  than  in  the  continental  mode  set.  This  does  not  give  complete 
coverage  in  the  period  range  considered,  but  should  give  an  adequate  Green’s  function. 

To  illustrate  the  effect  that  the  oceanic  propagation  has  on  the  FE  results  sam¬ 
pling  a  wavefield  that  has  already  propagated  through  a  transition  region  three  sample 
distances,  from  FE  node  to  receiver,  of  one  hundred,  two  hundred  and  fifty,  and  one 
thousand  kilometers  were  chosen.  Mode  sum  Green’s  functions  for  receivers  at  the 
crustal  surface  and  sources  where  the  FE  displacements  and  stresses  were  recorded 

were  calculated  to  be  used  in  the  RT  integral.  The  seismograms  resulting  from  the 

RT  integrations  are  illustrated  in  Figure  19.  In  this  figure  the  RT  integration  results 
for  the  two  shortest  paths  are  compared  to  FE  results  at  the  same  locations  with 
respect  to  the  source.  The  result  for  the  longer  path  is  also  shown  below  the  upper 
two  pairs  of  seismograms.  In  each  of  the  pairs  of  illustrated  seismograms  the  upper 
trace  shows  the  FE  result.  This  result  includes  propagation  through  an  oceanic  path 
length  of  125  km,  for  the  upper  pair,  or  275  km,  for  the  lower  pair.  These  seismo¬ 
grams  are  taken  from  the  reverse  reference  model  FE  calculation  which  uses  a  column 
of  nodes  recorded  in  the  fifty  kilometer  forward  transition  calculation  as  forcing  func¬ 
tions.  Displacements  and  stresses  recorded  at  the  same  horizontal  distance  from  the 

source  as  these  forcing  functions  are  used  as  the  u2,  and  mu2i,  terms  in  the  RT 
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integration.  The  seismograms  resulting  from  the  RT  integration  of  these  terms  with 
Green’s  functions  for  propagation  through  a  path  length  of  oceanic  structure  are 
shown  as  the  lower  seismograms.  The  upper  pair  of  seismograms  is  for  Green’s 
functions  for  a  one  hundred  kilometer  path  length,  the  lower  pair  for  a  Green’s 
function  path  length  of  two  hundred  fifty  kilometers.  The  numbers  to  the  right  of, 
each  pair  indicate  the  RMS  amplitude  ratio  of  the  upper  seismogram  to  the  lower  one. 
The  upper  pair  of  seismograms  shows  reasonable  agreement.  At  least  half  the 
difference  between  the  two  traces  could  easily  be  due  to  truncation  errors  caused  by 
terminating  the  integration  surface  at  35  km  depth.  The  remaining  differences  in 
amplitude  and  waveform  can  be  explained  as  being  due  to  a  small  component  of  non- 
modal  energy  which  is  filtered  out  of  the  result  when  using  the  RT  integration  but  not 
when  using  the  FE  method.  Thus,  this  pair  of  seismograms  shows  that  after 
propagation  through  125  km  of  oceanic  path  most  of  the  non-modal  energy  is  no 
longer  visible  at  the  surface.  From  the  FE  calculations  whose  time  slices  are  shown  in 
chapter  two  it  is  clear  that  at  this  distance  some  of  the  non-modal  energy  is  still 
present  at  depth.  This  energy  is  probably  a  source  of  the  amplitude  discrepancy 
between  the  two  uppermost  seismograms.  The  lower  pair  of  seismograms  shows 
excellent  agreement  in  amplitude  but  discrepancies  in  waveform  for  all  but  the  initial 
arrivals.  This  is  expected  since  the  FE  result  is  contaminated  with  small  reflections 
from  the  nearby  grid  edge  that  can  change  the  waveform  significantly  but  have  been 
shown  to  change  the  RMS  amplitudes  by  at  most  a  couple  of  percent.  Thus  the 
amplitude  agreement  at  this  distance  implies  that  the  non-modal  energy  remaining  at 
depths  included  in  the  integration  is  not  significant.  However,  the  remaining 
amplitude  of  the  Lg  phase  is  still  fifty  five  to  seventy  percent  of  the  incident 
amplitude.  Using  longer  transitions  this  can  be  reduced  another  five  to  ten  percent, 
and  adding  the  effects  of  intrinsic  or  scattering  attenuation,  which  are  not  considered 


V 
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in  the  calculations  presented  here,  might  produce  a  further  small  reduction  in  the  the 
Lg  wave  amplitude.  Thus,  it  is  apparent  that  the  attenuation  due  to  passage  through 
a  transition  region  is  not  sufficient  to  explain  the  entire  observed  attenuation  of  Lg 
waves. 

In  conclusion,  the  FE  results  show  several  effects  which  are  significant  mechan-- 
isms  for  the  attenuation  of  Lg  as  it  passes  through  ocean  continent  or  continent  ocean 
transition  regions.  The  most  important  of  these  effects,  associated  with  the  path 
length  of  oceanic  structure,  is  the  transmission  of  energy  across  the  crust  half-space 
boundary  in  the  forward  transition  region  and  the  conversion  of  that  energy  into 
downgoing  non-modal  phases.  As  the  length  of  the  oceanic  path  increases  the  ampli¬ 
tude  of  these  downgoing  phases  increase  and  they  are  able  to  propagate  far  enough 
towards  the  base  of  the  grid  that  they  are  not  transmitted  back  into  the  crustal  layer 
through  the  half-space  crust  boundary  in  the  reverse  transition  region.  The  FE  results 
illustrated  in  this  chapter  show  that  significant  portions  of  the  energy  originally 
trapped  in  the  continental  crust  can  be  carried  out  of  the  system  in  this  manner.  The 
RT  integration  results  show  that,  for  a  125  km  path  length,  much  of  the  downgoing 
energy  has  passed  out  of  the  depth  region  above  the  base  of  the  continental  crust,  and 
that  after  propagation  through  an  oceanic  path  of  length  275  km  the  non-modal 
energy  has  become  negligable.  The  attenuation  of  Lg  passing  through  the  reverse 
transition  should  be  maximal  for  path  lengths  in  excess  of  that  length  at  which  the 
non-modal  energy  has  reached  depths  below  the  base  of  the  continental  crust.  The 
oceanic  path  length  necessary  for  complete  escape  of  the  non-modal  component  of  the 
wavefield  is  on  the  order  of  200  to  250  km  but  significant  portions  of  the  energy  is 
able  to  escape  for  path  lengths  as  short  as  eighty  kilometers.  The  results  presented 
here  are  for  a  fifty  kilometer  long  transition.  Since  a  longer  transition  region  produces 
more  rapid  downward  propagation  these  lengths  should  be  expected  to  be  an 


overestimate  of  the  observed  Lg  data  which  is  oq  average  associated  with  longer  tran¬ 
sitions.  However,  the  amount  of  attenuation  seen  in  the  Lg  phase  is  not  large  enough 
to  explain  the  observed  attenuation  of  Lg.  These  simple  models  show  that  geometry 
of  the  transition  region  alone  produces  at  most  a  reduction  of  a  factor  of  five  or  six  in 
amplitude.  Thus,  the  observed  attenuation,  which  has  frequently  been  attributed  to 
the  propagation  effects  caused  by  passage  through  a  transition,  cannot  be  entirely  due 
to  these  propagation  effects. 
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